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FEASIBILTY  OF  A  BLAST  WAVE 
ATTENUATION  STRUCTURE 


by  Dale  Richard  Hartmann 

Chairman  of  the  Supervisory  Committee:  Professor  Ashley  F.  Emery 
Department  of  Mechanical  Engineering 


This  thesis  begins  with  an  overview  of  bombings  in  the  United  States,  followed  by  the 
introduction  of  the  Rankine-Hugoniot  equations  for  blast  wave  pressure.  The  subsequent 
chapters  develop  the  one-dimensional  and  two-dimensional  Euler  equations.  These 
equations  are  the  solved  using  the  MacCormack  finite  difference  algorithm.  The  basis  of 
the  investigation  then  begins  by  placing  pole,  shear  plate  and  wedge  obstacles  in  the  path 
of  the  blast  wave.  The  results  of  these  simulations  are  interpreted  and  conclusions 
presented.  Finally  a  synopsis  of  the  existing  results  and  cost  analysis  for  structure 
hardening  are  presented. 
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PREFACE 


The  inspiration  for  this  thesis  was  developed  over  the  course  of  two  years.  While  at  an 
assignment  with  the  Defense  Nuclear  Agency,  I  was  privileged  to  participate  in  several 
tests  involving  car  bombs.  I  became  very  interested  with  the  effects  of  the  blast  waves  on 
different  types  of  structures.  After  familiarizing  myself  with  the  current  methods  of 
protecting  an  existing  structure  from  blast  waves,  I  thought  there  must  be  a  better  way. 
This  thesis  is  my  investigation  regarding  whether  there  is  a  better  way  to  protect  existing 
structures. 
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INTRODUCTION 


BOMBS:  THE  UNKNOWN  MENACE 


The  threat  of  terrorist  bombings  is  present  every  year.  Although  the  United  States  has 
been  relatively  immune  from  terrorist  bombings  such  as  those  in  Ireland  or  Israel,  2577 
bombings  occurred  in  the  United  States  in  1995  alone.  Bombings  involving  improvised 
explosive  devices  such  as  pipe  bombs  numbered  1,562  in  1995.  The  majority  of  the 
bombings  in  1995  were  small  in  size  but  produced  $105  million  damage  and  937 
casualties  (193  killed  and  744  injured).  The  worst  bombing  in  the  history  of  the  United 
States,  the  Murrah  building  in  Oklahoma  City,  accounted  for  $100  million  damage  and 
786  casualties  (168  killed  and  518  injured). 

A  cursory  review  of  the  above  facts  and  events  leads  to  the  question:  How  can  buildings 
that  are  possible  terrorist  bombing  targets  be  protected  from  the  effects  of  bomb  blasts? 
That  question  is  the  basis  of  this  thesis. 


CHAPTER  1 :  THE  PROBLEM 


TYPES  OF  EXPLOSIVES 

Explosives  classified  as  either:  High  explosives  or  Low  explosives.  High  explosives 
possess  certain  characteristics  that  are  vastly  different  than  Low  explosives.  High 
explosives  detonate,  which  means  that  when  they  are  initiated  a  shock  or  blast  wave  will 
form.  They  also  will  burst  or  shatter  materials  near  them,  are  capable  of  penetrating 
materials,  and  have  the  capability  of  lifting  or  moving  objects. 

Low  explosives,  in  contrast,  do  not  detonate  but  bum  very  rapidly.  Since  Low  explosives 
do  not  detonate,  the  pressure  rise  that  is  produced  is  usually  smaller  in  amplitude  but 
longer  in  duration  than  that  of  a  High  explosive.  This  combination  tends  not  to  produce 
an  impulse  type  of  blast  wave  but  a  slightly  'softer'  shock  to  materials  nearby. 

BLAST  WAVE  CHARACTERISTICS 

A  blast  wave  generated  as  the  result  of  the  initiation  of  a  contained  High  explosive  is 
created  in  the  following  sequence  of  events.  First  hot  gases  with  temperatures  of  the  order 
3000  degrees  C  are  generated.  In  concert  with  this  temperature  rise  the  pressure  of  these 
gases  is  of  the  order  100  to  300  kilobars.  It  is  this  second  characteristic  of  high  explosive 
initiation  that  is  of  concern  here. 

This  hot  high-pressure  gas  then  expands  into  the  surrounding  atmosphere.  As  expansion 
occurs,  the  air  surrounding  the  expanding  gas  is  forced  outward.  Since  air  is  a 
compressible  medium,  a  layer  of  air  adjacent  to  the  outer  edge  of  the  expanding  gas  is 
compressed.  It  is  this  layer  of  compressed  air,  which  is  called  a  blast  wave.  As  the  hot 
gas  expands  and  cools,  its  pressure  falls.  The  pressure  of  the  layer  of  compressed  air  also 
falls  as  it  is  pushed  farther  outward  from  the  point  of  detonation.  As  the  gas  continues  to 


cool  and  expand,  at  some  point  in  time  and  space  the  pressure  will  fall  below  atmospheric 
pressure  due  to  the  momentum  of  the  layer  of  compressed  air.  This  slight  negative 
pressure  is  called  overexpansion  and  results  in  a  negative  phase  of  the  blast  wave 
whereby  the  outward  flow  of  gases  and  air  is  reversed.  Eventually  after  some  number  of 
oscillations,  equilibrium  is  reached  and  the  motion  of  the  air  stops. 

The  speed  at  which  the  blast  wave  travels  is  different  depending  on  the  medium  in  which 
it  is  traveling.  For  a  bomb  resting  on  the  ground,  two  waves  will  exist:  one  traveling 
through  the  ground  and  another  traveling  through  the  air.  The  latter  is  what  this  thesis  is 
concerned  with. 

The  magnitude  of  the  impulse  of  an  explosive  blast  wave  is  dependent  upon  many 
factors.  Some  of  these  factors  are: 

Is  the  explosive  cased? 

What  type  of  explosive  is  used? 

The  quantity  of  explosive  used? 

What  weather  conditions  are  present  at  the  time  of  detonation? 

What  is  the  terrain? 

Are  any  structures  nearby? 

PROTECTING  EXISTING  STRUCTURES 

THE  TRADITIONAL  APPROACH 

The  traditional  approach  to  protecting  an  existing  structure  from  explosive  blast  damage 
is  to  harden  it,  that  is  to  increase  the  structural  strength  of  single  or  multiple  components. 
Hardening  of  a  building  consists  of  several  methods. 


The  first  and  most  rudimentary  method  is  to  harden  the  glazing,  if  any  is  present,  on  the 
exterior  of  the  structure.  Even  a  relatively  small  bomb  can  cause  large  amounts  of  injury 
and  death  without  structurally  damaging  a  building.  This  carnage  is  achieved  by 
shattering  normal  window  glass  and  propelling  the  glass  fragments  at  high  speed 
throughout  the  interior  of  rooms  along  the  outside  perimeter  of  the  building,  anyone 
present  in  the  rooms  is  impaled  with  large  shards  of  glass.  The  use  of  tempered  or 
tempered  safety  glass  can  minimize  this  effect.  However,  the  use  of  blast  curtains  (heavy 
polymer  curtains  weighted  along  the  lower  edge)  or  the  installation  of  Mylar  film  with 
reinforced  window  frames  can  eliminate  this  hazard.  The  problem  is  the  curtains  have  to 
be  kept  closed  at  all  times  to  provide  protection,  thus  rendering  the  window  useless  as  an 
architectural  entity.  In  addition  the  Mylar  is  relatively  expensive. 

The  most  complex,  and  best,  method  is  to  harden  the  entire  structure  against  the  effects  of 
explosive  blast  waves.  This  generally  consists  of  a  major  reconstruction  of  the  interior 
structural  elements  of  a  building  or  encasing  the  building  inside  of  another  more  robust 
exterior  structure  capable  of  withstanding  the  blast  wave.  This  method,  although  highly 
effective,  is  also  difficult  to  design,  disruptive  to  the  occupants,  time  consuming  to 
construct  and  very  expensive. 

THE  NEW  APPROACH 

The  approach  that  this  thesis  investigates  is  the  use  of  a  blast  wave  attenuation  structure. 
The  form  of  the  blast  wave  attenuation  structure  would  be  such  that  as  the  blast  wave 
travels  through  it  the  energy  of  the  wave  is  dissipated.  Ideally,  the  energy  of  the  wave 
would  be  lowered  to  such  a  degree  as  to  produce  and  overpressure  of  no  more  than  one 
half  atmosphere.  This  lowered  pressure,  although  still  high  enough  to  cause  some 
damage,  is  safe  for  most  commercial  types  of  construction. 

The  forms  of  the  attenuation  structure  to  be  investigated  include  shear  plates,  a  field  of 
poles,  and  of  wedges.  The  success  of  a  type  of  structure  will  depend  upon  its  ability  to 
attenuate  the  pressure  in  less  than  15  meters,  constructability,  low  cost,  and  aesthetics. 


CHAPTER  2:  THE  QUESTION 


MODELING  BLAST  WAVES 

BLAST  WAVE  PARAMETERS 

The  modeling  of  a  blast  wave  must  address  two  issues.  The  first  is  the  modeling  of  the 
static  aspects  of  the  explosion.  These  static  aspects  include  pressure,  density, 
temperature  and  shock  wave  velocity.  The  second,  and  more  difficult,  is  the  modeling  of 
the  dynamics  of  the  blast  wave  motion  and  interactions  with  objects.  These  dynamic 
aspects  are  the  same  quantities  as  the  static  aspects  but  are  concerned  with  how  these 
quantities  change  with  time  and  position. 

STATIC  ASPECTS 

It  is  important  to  model  the  static  aspects  of  an  explosion  since  this  is  what  provides  the 
driving  force  for  the  dynamic  solutions  and  simulations.  The  first  static  aspect  to  be 
modeled  is  the  pressure  generated  by  the  detonation  of  the  explosive.  The  pressure 
generated  is  dependent  upon  several  factors: 

Type  of  explosive 

Distance  from  explosive 

Position  of  explosive  relative  to  the  ground 

The  type  of  explosive  is  important  since  each  different  type  of  explosive  contains  a 
different  amount  of  energy  per  unit  mass.  These  differences  are  summarized  in  the  table 
below. 


Table  1  Explosive  Equivalents 


Explosive 


RDX  (Cyclonite) 
HMX 

Nitroglycerin  (Liquid) 

TNT 

Blasting  Gelatin 

Nitroglycerin  dynamite 

Semtex 

Compound  B 


Mass  Specific  Energy 

TNT  Equivalent 

Qx(kJ/kg) 

(Qx/Qtnt) 

5360 

1.185 

5680 

1.256 

6700 

1.481 

4520 

1.000 

4520 

1.000 

2710 

0.600 

5660 

1.250 

5190 

1.148 

As  can  be  seen  from  the  table  above,  the  method  of  normalizing  the  energies  of  different 
explosives  relative  to  that  contained  in  TNT  has  been  developed  and  is  universally 
accepted. 

The  relationship  between  the  range  and  the  TNT  equivalent  charge  mass  is  known  as  the 
scaled  distance  and  is  given  by  the  following  equation: 


R 


z  =  ^u;  (i) 


Where  W,  the  universal  symbol  for  TNT  equivalents,  denotes  the  mass  of  explosive  and 
R  denotes  the  range  from  the  detonation  of  the  explosive  to  the  point  of  interest.  Rankine 
and  Hugoniot  produced  the  first  analytic  solutions  for  blast  wave  front  parameters  in 
1870  for  normal  shocks  in  ideal  gases.  These  solutions  were  later  expanded  by  Brode  to 
determine  the  peak  static  overpressure  in  the  near  and  medium  fields  of  the  blast  wave. 
The  near  field  is  the  region  where  the  overpressure  is  higher  than  10  bar  (10  Pa).  The 
medium  field  is  the  region  where  the  overpressure  is  less  than  10  bar  (106  Pa)  and  higher 
than  0. 1  bar  (10   Pa).  The  equation  for  the  overpressure  in  the  near  field  is  given  by 

A=(fr+1)  OOOOOO)  (Pa)  (2) 

The  equation  for  the  overpressure  in  the  medium  field  is  given  by 


(0.975      1.455     5.85  V  x 

p,  =^-— -  +  —3-  +  _-0.019j(l00000)  (Pa)  (3) 

The  absolute  pressure  of  the  detonation  is  given  by 

*    absolute  is         r ambient  \     / 

The  equations  developed  by  Rankine,  Hugoniot  and  Brode,  shown  above,  dealt  with 
shock  waves  in  free  air.  These  shock  waves  were  allowed  to  propagate  in  all  directions, 
or,  more  simply,  in  a  spherical  manner.  Since  I  am  concerned  with  terrorist  bombings, 
the  majority  of  which  involve  explosives  placed  on  or  near  the  ground,  a  correction  factor 
is  needed  on  pabsoiute  to  account  for  the  reflection  from  the  ground.  If  the  ground  were  an 
ideal  surface  for  reflection  this  correction  factor  would  by  2.  However,  the  ground  will 
respond  in  two  ways  that  are  not  ideal.  The  ground  underneath  the  explosion  will 
compress.  The  area  around  the  explosion  will  undergo  a  fracturing  process,  which  will 
result  in  particles  being  ejected  into  the  air.  There  is  no  analytical  method  of  quantifying 
these  responses,  but  empirical  evidence  suggests  a  correction  factor  given  by 

Pa    =    X&P  absolute  (5) 


This  correction  factor  is  applied  to  the  pressure  calculations  for  an  explosion  to  provide 
the  input  pressure  for  all  analytical  and  finite  difference  calculations.  The  Mach  number 
of  the  detonation  shock  wave  is  given  by 


ff  v  .  ,  1  \V2 


M.  = 


^--l 


V  V  P ambient  ' 


y  +  \ 

K2y  +  \) 


(6) 


The  units  in  the  equations  below  are  correct  for  pa  expressed  in  Pa. 
The  density  of  the  gas  behind  the  detonation  Shockwave  is  given  by 


Pi 


r  ambient 

1- 

V 

2    \ 

r  +  v 

f>  '  1 

I       M  2) 

(Kg/m3) 


(7) 


The  velocity  of  the  detonation  Shockwave  is  given  by 


f    2    /- 


«■>  =  c, 


A/ 


(m/s) 


(8) 


Where  ci  is  the  ambient  speed  of  sound 


r  ambient 


1/2 


Y 

^       r  ambient  ' 


(m/s) 


(9) 


The  temperature  of  the  gas  behind  the  detonation  Shockwave  is  given  by 


T-,  =  -^-  (K) 
"      p2R 


(10) 


Where 


(11) 


With 


Ru  =  8314.51(J/kmoleK) 


MWair  =  29  (Kg/kmole) 


(12) 
(13) 


The  speed  of  sound  behind  the  detonation  Shockwave  is  given  by 


C-,  -    c. 


T2 

rj-i 

ambient  ^ 


1/2 


(m/s) 


(14) 


DYNAMIC  ASPECTS 


THE  EULER  EQUATIONS 

Air  has  a  small  thermal  conductivity  (k)  and  also  a  small  viscosity  (o)  therefore  the  terms 
in  the  Navier-Stokes  equations  that  depend  on  k  and  o  can  be  neglected,  which  leads  to 
equations  in  which  the  convective  terms  dominate  and  the  fluid  (in  this  case  air)  is  treated 
as  inviscid.  The  inviscid  2-D  Navier-Stokes  equations  are  called  the  Euler  equations  and 
are  the  equations  dealt  with  in  this  paper.  These  equations  are  much  simpler  in  form  and 
also  programming  effort. 

The  Euler  equations  are  written  as 


dQ     cF     dG 

dt       dx.      cty 


(15) 


where 


0  = 


P 

pu 
pv 
E 


F  = 


pu 
pu2  +  p 

puv 
u(E  +  p) 


G  = 


pv 

puv 
pv2  +p 
v(E  +  p) 


(16) 


and 
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^  +  I(M2+v2)J  (17) 


The  Euler  equations  represent  a  system  of  equations  that  conserve  mass,  momentum  and 
energy.  These  equations  must  be  solved  as  a  system  to  capture  nonlinearities  such  as 
shocks. 

The  Euler  equations  also  need  an  equation  of  state,  which  is  the  ideal  gas  law. 

p  =  pRT  =  (r-l)[E-^p{u2+v2)j  (18) 


CHAPTER  3:  ONE  DIMENSIONAL  EQUATIONS 

THE  1-D  EULER  EQUATIONS 

The  first  step  in  modeling  the  blast  wave  is  to  construct  a  one-dimensional  model  with 
time  and  distance  as  the  parameters.  The  1-D  Euler  equations  are  written  as 


—  + —  =  0 
dt      3c 


(19) 


where 


Q 


p 

pu 
E 


F  = 


pu 
pu"  +  p 
u(E  +  p) 


(20) 


and 


E  =  p[e  +  -{u2) 


(21) 


The  equation  of  state  is 


p  =  pRT  =  (r-l)  \e--p[u2) 


(22) 


This  system  of  equations  yields  the  following  three  equations 


Conservation  of  mass 


dp  |  o(pu)  =  Q 
d        3c 


(23) 
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Conservation  of  momentum 


oipu)     4pu2  +  p) 
yrj+^ ;-  =  0  (24) 


dt  3c 


Conservation  of  energy 


cE      ME  +  p) 

—  +  — * ^-  =  0  (25) 

dt  3c 


FINITE  DIFFERENCING 

To  evaluate  this  system  of  equations  the  method  of  finite  differencing  is  used.  Finite 
differencing  is  a  powerful  numerical  technique  that  can  be  used  to  solve  partial 
differential  equations.  Finite  differencing  algorithms  come  in  two  varieties:  explicit  and 
implicit.  An  explicit  algorithm  solves  each  equation  in  turn  and  then  applies  those  results 
to  the  initial  values  for  the  next  time  step.  An  implicit  algorithm  solves  the  system  of 
equations  simultaneously  via  matrix  manipulations.  I  chose  the  explicit  form  of  finite 
differencing  instead  of  the  implicit  form  for  the  following  reasons: 

Ease  of  coding 

Lower  computational  requirements 

Ability  to  run  on  desktop  computing  systems 

The  first  step  in  constructing  a  finite  differencing  program  is  to  choose  an  algorithm. 
Many  different  algorithms  exist  but  to  be  useful  for  the  problem  at  hand  the  chosen 
algorithm  must  be  able  to  "capture"  the  shock  discontinuities  well  and  also  handle 
surface  interactions.  Several  algorithms  meet  these  requirements: 

Forward  Time  Centered  Space  (FTCS) 

MacCormack 

Harten-Ye 
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1  st  Order  Roe 

I  initially  investigated  the  use  of  the  FTCS  algorithm,  however  I  was  not  pleased  with  the 
shock  "capturing"  abilities  of  this  algorithm.  This  algorithm  smears  or  spreads  the  shock 
discontinuity  over  4  to  5  Ax  intervals,  it  also  has  dispersion  errors  which  appear  as 
oscillations  at  the  shock  discontinuity.    The  following  figure  demonstrates  these  points. 
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Figure  1  FTCS  Output 


The  next  algorithm  I  investigated  was  MacCormack's.  This  algorithm  is  a  two-step 
method  that  is  known  for  "capturing"  shock  discontinuities  very  well.  In  contrast  to  the 
FTCS  algorithm,  MacCormack's  does  not  have  the  dispersion  errors  at  the  shock 
discontinuity.  The  following  figure  demonstrates  these  points. 
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Figure  2  MacCormack  Output 


The  Harten-Yee  and  1st  Order  Roe  methods  were  abandoned  early  on  due  to  their 
complexity  and  coding  difficulties  and  because  they  have  been  shown  not  to  be 
substantially  superior  to  MacCormack's  method. 

Due  to  the  simplicity  and  shock  capturing  abilities  of  the  MacCormack  algorithm,  I  have 
chosen  to  use  it  as  the  basis  for  the  remainder  of  this  thesis. 

The  method  of  the  MacCormack  algorithm  is  contained  in  the  two  steps  below: 


Predictor  step:  Qj  =  Q"  -  AtAxF" 


1/—     „„\     At„  — 


Corrector  step:  Q.n+]  =  -{q.  +  Q"\ V^  F. 


(26) 
(27) 


where 
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K  =  f.+;-f; 


(28) 


V    =  F  -  F 


(29) 


For  the  one  dimensional  case  the  MacCormack  algorithm  equations  become 


Conservation  of  mass 


dp  [  c(pu)  =  Q 
d        dx 


(30) 


Predictor:  Pj  =  p.n    -  -j-  i^pu  )  j+"  -  (pu)."j 


(31) 


Pj-Pj        M  l- 


CoiTector:  p."     =  — \pu  "  -  pu;  , " 

2  2Ajc\     ;  ' 


(32) 


Conservation  of  momentum 


c(pu)     o{pir  +p) 
dt  dx 


(33) 


Predictor:  (pu) .  =  (pu)  " 


Ax 


(H,,,j2  H")2 


n  n 


Pj< 


Pj 


+  Pj+X      ~Pj 


(34) 


i      (/-**)  ■  ~~  0°")  •"       A?  f  / — \      / — \         —      —    > 
Corrector:  (pu)"+   = J— '—  -  ^{[Pu).  "  (P") ._,  +  P;  -  Pj-i)       (35) 


Conservation  of  energy 


Predictor:  E .  -  E" 


Ax 
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cE     di(E  +  p) 

—  +  — ^ ^  =  0  (36) 

dl  dx 


\  r_j+\ 


E .  —  E "       At  i —      —         —      —     \ 
Corrector:  £."+1  =  — ^ J- [Ej  -  Ej-i+p,-p,  .)  (38) 

2  2Ajc^  j        J    ' 

The  equations  for  the  conservation  of  momentum  and  energy  above  contain  a  pressure 
term  that  is  evaluated  using  the  equation  of  state.  The  equation  of  state  is  the  ideal  gas 
law  and  is  given  by 


p  =  (r-i) 


_'W 


2     p 


(39) 


In  the  above  form  the  equation  state  can  be  evaluated  using  the  results  of  the 
MacCormack  algorithm. 

The  final  concern  for  any  finite  differencing  algorithm  is  stability.  The  use  of  the 
Courant-Friedrichs-Lewy  or  CFL  number  is  the  most  widely  used  method  of  controlling 
stability.  The  CFL  number  arises  out  of  the  results  from  a  Von-Neumann  stability 
analysis.  The  theory  of  Von-Neumann  utilizes  a  Fourier  transform  to  transform  the  finite 
difference  solution  into  a  wave  space  solution.  The  amplitudes  of  the  waves  in  this 
resulting  solution  will  grow  or  decay  based  upon  the  particular  finite  difference  algorithm 
chosen  and  the  value  of  the  CFL  number.  The  waves  in  the  Von-Neumann  solution  that 
grow  are  unstable,  those  that  decay  are  stable.  The  results  of  the  Von-Neumann  analysis 
require  that  the  CFL  number  be  less  than  or  equal  to  one.  The  physical  implication  of  the 
CFL  number  being  less  than  or  equal  to  one  is  that  the  sum  of  the  amplitudes  of  the 
waves  that  decay  is  larger  than  the  sum  of  the  amplitude  of  the  waves  that  grow.    Thus  if 
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the  Courant  number  remains  below  the  value  of  1  the  decaying  waves  will  dominate  and 
the  algorithm  will  remain  stable. 

Ax 

CFL  =  — r—  (40) 

\\u\  +  c)At 

where  c  is  the  speed  of  sound  for  the  given  pressure  and  density  and  CFL  is  the  Courant- 
Friedrichs-Lewy  number.  To  determine  the  time  step  used  in  the  finite  differencing  code 
the  CFL  equation  above  is  manipulated  such  that  the  following  equation  for  a  one- 
dimensional  algorithm  time  step  results. 

At  =  CFL-^  (41) 


VALIDATING  THE  ALGORITHM 

At  this  point,  after  developing  the  equations  for  the  one-dimensional  case  and 
constructing  the  finite  difference  algorithm,  validation  of  the  method  is  appropriate.  To 
validate  the  code  a  model  was  constructed  with  a  right  traveling  shock  that  reflects  from 
an  infinite  wall  at  the  right  most-boundary.  I  selected  this  validation  method  due  to  its 
ease  of  coding  and  readily  available  analytic  solution. 
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Figure  3  Shocktube  with  right  traveling  wave 


The  actual  validation  was  done  by  allowing  the  finite  difference  solution  to  run  long 
enough  following  the  interaction  with  the  right  most  wall  and  then  comparing  the 
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pressure  and  density  values  from  the  model  to  the  analytic  values  calculated  from  the 
analytic  solution  (Shapiro) 

The  constant  initial  conditions  in  the  program  code  are 

Table  2  Fixed  initial  conditions 


Temperature        Gamma  Density  Pressure         Speed  of  sound 

298  K  i~4  1.614  kg/mJ       101325  Pa  296  m/s 


The  initial  conditions  input  to  the  finite  difference  model  were 

Table  3  Input  initial  conditions 


Courant 

Calculation 

Number  of 

Number  of 

Distance 

Explosive 

Explosive 

number 

distance 

position 
steps 

time  steps 

from  charge 

mass 

0.5 

4  m 

101 

250 

10m 

TNT 

200  kg 

The  program  calculated  values  behind  the  shock  prior  to  the  wall  interface  were 
determined  using  the  Rankine-Hugoniot  and  Brode  equations,  equations  2,  3,  5,  6,  7,  8,  9, 
10  and  14.  The  values  output  from  the  program  are  tabulated  below. 
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Table  4  Calculated  free  stream  shock  values 


Static  pressure  of  blast 

5.7417x10"  Pa 

Mach  number  of  shock 

2.2504 

Density 

4.8729  kg/mJ 

Shock  velocity 

446.18  m/s 

Temperature 

410.97  K 

Speed  of  sound 

348.97  m/s 

The  calculated  values  from  the  program  behind  the  shock  after  the  wall  interface  were 


Table  5  Calculated  reflection  shock  values 


Static  pressure  of  blast         2.246x10  Pa 


12 .17 '9  kg/mJ 


Density 


To  obtain  the  analytic  values  an  iterative  method  was  used.  Following  the  reflection 
from  the  right  most  wall  the  shock  is  now  traveling  to  the  left.  The  analytic  solution  is  to 
hold  the  shock  stationary  and  iterate  for  the  value  of  W. 
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Figure  4  Shocktube  with  stationary  wave 


The  equations  required  for  the  iterative  solutions  are 


M 


w  +  u. 


(42) 


£-=1 

Pi 


2Y 

y  + 1 


M 


(43) 


Ei. 

P2 


1- 


x  +  1 


M  "J 


(44) 


f    2    ' 


IA1-1  Li/  *} 


lr  +  i 


M  - 


MJ 


(45) 


Using  the  equations  above,  equations  41  through  44,  with  the  initial  values  listed  in  the 
table  below,  the  results  of  the  analytic  solution  will  be  those  listed  in  the  table  below. 
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Table  6  Analytic  results  for  shocktube 


Initial  values  u;  C2  p:  p2  y 

446.18  m/s  348.1525  m/s       574.17  kPa  4.8729  kg/m3         1.4 

Calculated  values      W  P3  p3 

180.5  m/s  2.0747  MPa  1 1.4963  kg/m3 


Comparing  these  analytic  values  with  the  finite  difference  calculated  values  yields  an 
error  of  +6%  for  the  density  and  +8%  for  the  pressure.  Which  means  that  the  finite 
difference  program  will  yield  a  slightly  conservative  (higher  than  true  values)  result. 
This  is  acceptable  for  my  purposes. 


CHAPTER  4:  TWO-DIMENSIONAL  EQUATIONS 


THE  2-D  EULER  EQUATIONS 

The  addition  of  a  second  dimension  to  the  1-D  Euler  equations  is  fairly  straightforward. 
The  second  dimension  requires  one  additional  equation  to  the  matrix  and  two  additional 
terms.  The  2-D  Euler  equations  are  written 


dQ     dF     SG     n 

—  +  —  +  —  =  0 

dt      3c      dp 


(46) 


where 


p 

pu 

pv 

pu 

F  = 

pu2  +  p 

G  = 

puv 

pv 

piiv 

pv2+p 

E 

u(E  +  p) 

_v{E  +  p) 

(47) 


and 


E  =  p[e  +  -[u2  +v2) 


(48) 


The  equation  of  state  is 


p  =  pRT  =  (r-\)^E~^p{u2+v2)j 


(49) 


This  system  of  equations  yields  the  following  four  equations 
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Conservation  of  mass 


*+fH  +  4H  =  0  ,50) 

d        3c  dy 


Conservation  of  x  momentum 


dt  3c  3y 

Conservation  of  y  momentum 


o(pu)  |  <{pu2+p)  |  <puv)  =  Q 


<fr») ,  4/H ,  4^2+^)  =  0  (52) 

<5  <^:  ^ 


Conservation  of  energy 


—  +  — * ^-  +  — * ^  =  0  (53) 

The  method  of  the  MacCormack  algorithm  for  two  dimensions  is  similar  to  that  for  one 
dimension.  The  two  dimensional  algorithm  is  contained  in  the  two  steps  below 

Predictor  step:  Q~k  =  QJk"  -  At(AxF,k"  +  A/7,/)  (54) 

Corrector  step:  g./+1  =  ±(Q~k  +  Q./)  -  ^(v J~k  +  Vy  G~7)  (55) 

where 

A,  =  Fj+X;  -  E/  (56) 


^x=F,k-F.hk  (57) 
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and 


a,  =  gjm;  -  gs; 


(58) 


Vy  =  Gjik  -  Ghk_x 


(59) 


For  the  two  dimensional  case  the  MacCormack  algorithm  equations  become 


Conservation  of  mass 


dp  |  d{pu)  |  c(pv)  =  Q 
d        3c  dy 


(60) 


Predictor:  pjk  =  pJk"  -  At 


(pu)j+ik"-{pu)jk"    {pv)jk+';-(pv). 


n\ 


Ax 


Ay 


(61) 


Corrector:  p . ; 


n  +  l 


£>;,*   -  /?;,*  A/ 


P",/  -  PUj-i/  +  {*j/  ZfhtL        (62) 


Ax 


Ay 


Conservation  of  x  momentum 


d{pu)     &(pu2  +  p)     d{puv) 
dt  dx  dy 


(63) 
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H„* =  H/ 


Predictor:  -A^ 


Ax 


(M,+u")      ((^)m") 


h  n 


Pj+Uk 


Pj.k 


+  /W     ~Pj,k 


(64) 


+  - 


Ay 


HwH^i"  W;/Wy, 


«A 


PjM\ 


Pj. 


(/*),, 


„+I  H,*-(H/ 


Corrector: 


A; 


Ax- 


HJ  -IH,--J  +^-/v>, 


(65) 


+■ 


Ay 


HJ^j,  H^H 


,*-i 


^y.* 


P 


/,*-! 


where 


P  =  (r-1) 


2     P 


2\ 


(66) 


and 


p  =  (r-0 


W 


2     p 


(67) 
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Conservation  of  y  momentum 


d(pu)     c{puv)     &(pv2  +  p) 


+  — -  + 


dt  3c 


$ 


(68) 


HrH/ 


Predictor:  -A? 


Ax 


„H;,i+l"W;,W"        H/  Wy, 


«A 


P** 


i,k  +  l 


P 


M  J 


(69) 


+■ 


Ay 


H+1/)2  H/f 


p> 


7  +  1,* 


P 


«  rJ 


+  Pm/-Pj/ 


J* 


H,, 


„,  H,,-(H,/ 


Corrector: 

2 


1 
Ajc 


H,-.*H,  H-uH-., 


Pyji 


P 


/-!,* 


(70) 


Ay 


+^IH,J  -IH.J  +^M-P^i 


where 
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and 


Conservation  of  energy 


p=(t-i) 


E     "W 


V 


2     p 


(71) 


p  =  (r-0 


£- 


R 


-\ 


/> 


(72) 


^     di(E+p)     d>(E+p)_ 


dt  3c 


& 


(73) 


^=V 


Predictor:  At 


Ax 


H/+u" 


H/ 


/>y+u 


—\&j+i,k     ~Pj+\,k   )  V ^y,*     ~".Pyf*   j 


(74) 


Py. 


+ 


Ay 


PjM\  Pj,k 
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,      E.-E" 


Coirector: 

2 


Ax 


v    //,*  /t;         -     \     v//-i,*  /-p  \ 

"j,k  Pj-l,k 


(75) 


f  I — 


+  - 


Ay 


R~„ 


"(^  +P;,,)-— ^±(^-1  +^-.) 


/> 


M-i 


The  pressure  terms  in  the  conservation  of  energy  equations  above  are  evaluated  by  using 
the  equation  of  state.  The  pressure  terms  are  distinct  for  each  direction.  The  x  direction 
pressure  term  is  given  by 


(  f  (       \2> 


p  =  (y + i) 


p 


(76) 


Similarly  the  y  direction  pressure  term  is  given  by 


p  =  (r  +  i) 


.  p  J) 


(77) 


The  time  step  in  two  dimensions  is  a  more  complicated  expression  given  by 


At  =  CFL 


+ 


1 
+  c\ — -  +  ■ 


Ax     Ay       V  Ax"      Ay' 


(78) 
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where  c  is  the  speed  of  sound  for  the  given  density  and  pressure  and  CFL  is  the  Courant- 
Friedrichs-Lewy  number. 

VALIDATING  THE  ALGORITHM 

Applying  reasoning  similar  to  the  one-dimensional  case,  an  input  pressure  was  applied  to 
the  x  edge  of  the  computational  domain.  The  effect  of  this  input  pressure  was  the  same 
as  the  calculations  preformed  by  the  one-dimensional  code.  Inputting  the  same 
parameters  into  both  the  2-D  and  1-D  programs  allowed  comparison  of  the  output 
matrices,  which  were  identical.  Therefore  the  2-D  code  is  validated  since  the  1-D  code 
has  been  proved  correct  and  acceptable.  Following  validation  of  the  2-D  code  for  a  linear 
pressure  wave  traveling  in  the  y  direction,  the  same  procedure  was  applied  for  a  wave 
traveling  in  the  x  direction.  The  same  results  were  obtained. 

The  results  of  this  validation  are  shown  in  the  tables  below 
Table  7  1-D  program  calculated  values 


Density 

4.8729 

3.0026 

1.7713 

1.6140 

1.6140 

Momentum 

2174.2 

974.5 

79.3 

0 

0 

Energy 

1.9395e6 

1.0087e6 

0.3144e6 

0.2533e6 

0.2533e6 

The  results  for  a  2-D  wave  traveling  in  the  x  direction  are  given  in  the  tables  below 
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Table  8  2-D  x-traveling  wave  program  calculated 
values 


Density 

4.8729 

3.0026 

1.7713 

1.6140 

1.6140 

4.8729 

3.0026 

1.7713 

1.6140 

1.6140 

4.8729 

3.0026 

1.7713 

1.6140 

1.6140 

4.8729 

3.0026 

1.7713 

1.6140 

1.6140 

4.8729 

3.0026 

1.7713 

1.6140 

1.6140 

Momentum 

2174.2 

974.5 

79.3 

0 

0 

2174.2 

974.5 

79.3 

0 

0 

2174.2 

974.5 

79.3 

0 

0 

2174.2 

974.5 

79.3 

0 

0 

2174.2 

974.5 

79.3 

0 

0 

Energy 

1.9395e6 

1.0067e6 

0.3144e6 

0.2533e6 

0.2533e6 

1.9395e6 

1.0067e6 

0.3144e6 

0.2533e6 

0.2533e6 

1.9395e6 

1.0067e6 

0.3144e6 

0.2533e6 

0.2533e6 

1.9395e6 

1.0067e6 

0.3144e6 

0.2533e6 

0.2533e6 

1.9395e6 

1.0067e6 

0.3144e6 

0.2533e6 

0.2533e6 

The  results  for  a  y  traveling  wave  are  in  the  tables  below 
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Table  9  2-D  y-traveling  wave  program  calculated 
values 


Density 

4.8729 

4.8729 

4.8729 

4.8729 

4.8729 

3.0026 

3.0026 

3.0026 

3.0026 

3.0026 

1.7713 

1.7713 

1.7713 

1.7713 

1.7713 

1.6140 

1.6140 

1.6140 

1.6140 

1.6140 

1.6140 

1.6140 

1.6140 

1.6140 

1.6140 

Momentum 

2174.2 

2174.2 

2174.2 

2174.2 

2174.2 

974.5 

974.5 

974.5 

974.5 

974.5 

79.3 

79.3 

79.3 

79.3 

79.3 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

Energy 

1.9395e6 

1.9395e6 

1.9395e6 

1.9395e6 

1.9395e6 

1.0067e6 

1.0067e6 

1.0067e6 

1.0067e6 

1.0067e6 

0.3144e6 

0.3144e6 

0.3144e6 

0.3144e6 

0.3144e6 

0.2533e6 

0.2533e6 

0.2533e6 

0.2533e6 

0.2533e6 

0.2533e6 

0.2533e6 

0.2533e6 

0.2533e6 

0.2533e6 

The  results  above  were  achieved  with  the  following  input  conditions: 


Table  10  2-D  validation  initial  conditions 


Time  step 

x-grid 

Range             to 
Charge 

Explosive  Type 

Explosive 
Mass 

5e-4  s 

5 

10m 

TNT 

200  kg 

CHAPTER  5:  TESTING  THE  NEW  APPROACH 


The  validation  of  the  2-D  code,  presented  in  the  previous  chapter,  was  performed  with 
lineai-  wavefronts  traveling  in  only  one  direction.  While  this  situation  might  occur  with  a 
very  large  explosion  at  a  great  distance,  in  general  this  is  an  unrealistic  case.  Therefore, 
the  use  of  symmetry  was  applied.  The  final  step  in  the  code  development  process  was  to 
place  the  explosion  in  one  comer  of  the  computational  realm  and  allow  the  wave  to 
expand  in  true  2-D  fashion.  Placement  of  the  explosive  in  this  corner  allows  for  two 
planes  of  symmetry  and  reduces  the  number  of  calculations  per  time  step  by  a  factor  of 
two.  This  reduction  of  calculations  is  also  beneficial  since  it  ignores  the  blast  wave  that 
travels  away  from  the  obstacles. 

The  initial  parameters  for  an  explosion  were  placed  into  one  corner  of  the  computational 
realm  and  allowed  to  expand  freely  with  no  obstacles  present.  The  results  obtained  were 
consistent  with  the  expected  outcome.  This  step  was  necessary  to  provide  a  test  of  the  2- 
D  code  in  the  true  two  dimensional  environment. 

OBSTACLE  BOUNDARY  CONDITIONS 

To  obtain  meaningful  results  with  obstacles  in  the  path  of  the  wave,  a  thorough 
understanding  of  the  correct  boundary  conditions  for  the  obstacles  is  needed.  The 
minimum  size  of  an  obstacle  is  limited  to  a  grid  of  3  by  3  computational  points.  This 
minimum  size  is  necessary  to  allow  for  the  specification  of  conditions  inside  the  obstacle. 
There  is  no  limit  on  the  maximum  obstacle  size,  although  it  obviously  cannot  exceed  the 
size  of  the  computational  realm.  The  figure  below  illustrates  a  minimum  size  obstacle. 
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Figure  5  Obstacle  boundary  conditions 


As  can  be  seen  in  the  figure  above  a  3  x  3  grid  is  the  smallest  obstacle  possible  to  provide 
for  the  specification  of  quantities  inside  an  obstacle.  The  figure  above  also  shows  what 
quantities  exist  at  the  points  on  the  boundaries.  The  quantities  that  exist  at  the  interior 
point  are  density  and  energy. 

For  inviscid  flow  on  rigid  surfaces  four  distinct  boundary  conditions  exist.  The  first  is 
that  the  component  of  the  x-momentum  normal  to  the  surface  is  zero.  Similarly,  the 
component  of  the  y-momentum  normal  to  the  surface  is  also  zero.  The  third  condition  is 
that  the  pressure  gradient  is  zero. 
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The  last  boundary  condition  is  that  velocity  component  normal  to  the  obstacle  surface  is 
zero. 

Another  assumption  needed  to  evaluate  all  parameters  is  constant  energy  inside  the 
obstacle,  E.  The  equation  for  energy  is: 

£=p^  +  I(M2+v2)J  (79) 

Since  at  this  point  the  velocities  are  known,  the  energy  and  the  pressure  can  be 
calculated.  This  equation  along  with  the  velocities  allows  for  calculating  p  on  the 
surface. 

The  pressure  at  the  surface  of  the  body  obtained  by  solving  the  ideal  gas  law  equation, 
which  is: 


p  =  pRT  =  (y-l)\E-±p(u2+v2)j 


(80) 
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RESULTS  OF  VARIOUS  ATTENUATION  STRUCTURES 

The  following  sections  detail  the  results  of  running  the  simulation  with  poles,  shear  plates 
and  wedges  for  obstacles. 

Before  reviewing  the  results  of  the  simulations  it  is  helpful  to  examine  the  flow  field 
results  without  any  obstacles  present.  All  of  the  results  to  be  presented  are  at  30  time 
steps.  The  program  does  utilize  variable  time  stepping  so  the  results  are  not  at  the  same 
total  time.  These  results  are  shown  in  the  figure  below. 


Absolute  Pressure 
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Figure  6  2-D  flow  field  results  without  obstacles 

The  figure  above  illustrates  two  points  unique  to  numerical  solutions.  The  first  point  is 
the  gradient  at  the  shockfiront.  The  gradient  of  the  shock  is  shown  as  the  lines  spaced 
closely  together.    These  lines  follow  the  curved  path  of  the  Shockwave  and  travel  at  the 
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speed  of  the  shock.  The  other  point  that  is  well  illustrated,  are  the  numerical  oscillations 
inherent  in  this  type  of  solution.  The  regions  outlined  behind  the  shockfront  are,  in  fact, 
numerical  oscillations  of  the  solution  and  do  not  exist  in  the  real  explosion 


POLES 

The  results  for  the  field  of  poles  were  disappointing  since  it  was  the  conceptual  basis  for 
this  thesis.  The  pole  obstacles,  in  the  figure  below,  were  set  at  3  by  3  grid  points  with 
infinite  height.  The  actual  physical  size  of  the  obstacles  varied  with  the  total  number  of 
grid  points  and  the  calculation  distance.  For  example  a  10  by  10  grid,  with  a  4  meter 
calculation  distance  generates  a  pole  that  is  1.2  maters  on  a  side.  It  is  very  important  to 
ensure  that  the  obstacle  size  is  reasonable.  In  the  above  example  to  make  the  pole 
dimensions  smaller  a  greater  total  number  of  grid  points  would  be  needed  or  a  smaller 
calculation  distance. 

Following  the  completion  of  many  simulations  the  effect  of  poles  on  reducing  the 
pressure  of  a  blast  wave  it  was  found  that  the  effect  is  minimal.  The  blast  wave  does 
compress  the  air  on  the  blast  side  of  the  pole  as  I  expected,  but  the  shadow  of  the  pole 
does  not  extend  far  enough  in  space  to  actually  cause  a  decrease  in  the  pressure  of  the 
wave  after  it  has  passed  through  the  entire  field.  The  cause  of  this  phenomenon  is  the 
low  viscosity  of  air.  The  calculations  that  are  used  for  the  simulation  are  inviscid,  which 
allows  the  air  to  flow  easily  around  a  small  object  such  as  a  pole.  For  the  effect  of  an 
obstacle  to  be  maintained  farther  downrange,  the  obstacle  must  be  of  a  size  large  enough 
to  produce  a  significant  wake  or  a  field  of  small  poles  set  in  a  grid  or  staggered  grid 
pattern  that  are  close  together.  The  field  of  poles  would  need  to  be  set  on  the  order  of 
their  diameter  to  be  effective.  Due  to  simulations  running  on  a  desktop  computer  system 
I  could  not  generate  a  field  of  poles  of  sufficient  density  to  demonstrate  this  directly.  I 
am  extrapolating  a  solution  based  upon  the  limited  observations  and  conclusions  I  made 
following  the  simulations.  The  following  figure  illustrates  this  result. 


37 


Absolute  Prvuurv 


M> 


Figure  7  Flow  field  results  for  pole  obstacles 

This  figure  shows  that  the  first  pole,  located  at  grid  point  (3,3)  has  a  very  large  pressure 
gradient  on  its  face.  This  gradient  is  produce  by  the  compression  of  the  air  between  the 
immovable  pole  and  the  moving  air  behind  This  compressed  air  contains  a  great  deal  of 
energy  that  remains  stagnant  upon  the  face  of  this  pole. 

The  figure  also  shows  that  following  passage  of  the  shock  through  the  pole  field  the 
gradient  has  been  spread  over  a  larger  area,  since  the  contour  lines  are  farther  apart,  but 
the  magnitude  of  the  gradient  has  not  changed.  Therefore,  the  shock  is  not  as  steep  as  it 
was  originally  but  is  of  essentially  the  same  magnitude.  Thus,  the  field  of  poles  has 
'softened'  the  shock  but  has  not  attenuated  the  peak  pressure.  The  interpretation  of  this 
result  is  that  this  spacing  of  poles,  which  is  an  example  of  many  simulation  runs,  does  not 
achieve  the  desired  outcome.  Extending  this  interpretation  further  leads  to  the  conclusion 
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that  no  economical  pole  field  would  attenuate  a  shock  to  a  level  low  enough  to  protect  a 
structure. 

SHEAR  PLATES 

A  shear  plate  is  a  thin  long  plate  with  the  long  axis  ideally  placed  parallel  to  the  direction 
of  flow.  Due  to  the  uncertainty  in  the  location  of  an  explosive  device,  I  placed  the  shear 
plates  such  that  the  length  axis  is  perpendicular  to  the  face  of  the  protected  building.  The 
shear  plates,  in  the  figure  below,  are  3  grid  points  wide  by  5  grid  points  long.  I 
performed  numerous  simulations  with  other  sizes  and  orientations  but  this  combination  of 
width,  length  and  placement  best  demonstrates  the  effect  of  shear  plates  on  the  shock. 

The  shear  plates  performed  slightly  better  than  the  pole  field  largely  due  to  their  size  and 
orientation.  The  first  plate  redirects  the  majority  of  the  flow  field  in  such  a  way  that  it  no 
longer  raises  the  pressure  on  the  lower  surfaces  of  subsequent  obstacles.  The  pressure  on 
the  upper  side  of  each  shear  plate  is  significantly  less  than  the  pressure  on  the  lower  side. 
This  is  the  result  of  the  lower  plates  acting  as  walls  and  shielding  the  upper  plates  from 
the  blast  pressure.  The  plates  also  redirect  the  flow  by  straightening  it  in  the  length 
direction  of  the  plates.  After  the  flow  has  passed  the  end  of  each  plate  it  begins  to 
diffract.  The  placement  of  the  shear  plates  is  critical  to  achieve  this  effect. 

The  results  of  the  simulation  can  be  seen  in  the  following  figure. 


39 


Ausolute  Pressure 


Figure    8    Flow    field    results    for    shear    plate 
obstacles 

The  figure  above  demonstrates  the  effect  of  shear  plates  on  the  flow.  The  lower  plate 
redirects  the  majority  of  the  flow  along  the  lower  edge  of  the  computational  realm  and 
removes  a  great  deal  of  the  energy  from  the  flow  that  impacts  the  middle  and  top  plates. 
The  flow  then  diffracts  around  the  lower  plate  after  it  has  passed.  The  gradient  is  again 
'softened',  but  its  magnitude  remains  essentially  the  same  as  the  original  unimpeded 
shock.  The  large  gradient  in  the  lower  left  corner  is  produced  by  the  shock  impacting 
upon  the  face  of  the  lower  shear  plate,  which  is  immovable.  This  impact  compresses  the 
air  and  produces  a  large  pressure  peak  in  a  very  small  area.  The  middle  and  top  shear 
plates  have  similar  gradients  on  their  faces  but  the  magnitude  is  much  smaller. 


To  provide  significant  protection  for  a  building  the  location  of  the  explosive  would  either 
need  to  be  known  or  surmised  before  construction  of  the  plates.  The  plates  could  then  be 
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placed  in  such  a  way  as  to  direct  the  flow  away  from  the  protected  structure. 
Unfortunately,  foresight  is  rarely  this  accurate  and  the  construction  of  many  plates  (or 
walls  for  larger  buildings)  will  be  neither  aesthetically  pleasing  nor  particularly 
inexpensive. 

WEDGES 

A  wedge  is  to  two  thin  plates  placed  at  an  angle  to  each  other  and  joined  at  the  apex.  The 
wedges  used  in  this  simulation  consist  of  two  plates  3  grid  points  by  5  grid  points  placed 
perpendicular  to  each  other.  The  wedges  were  then  placed  into  the  computational  realm 
such  that  the  lower  face  of  the  wedge  is  diagonally  across  the  path  of  the  shock.  The 
wedge  in  this  position  acts  somewhat  as  a  shield  or  wall  in  the  path  of  the  shock. 

The  results  for  the  wedge  obstacles  were  similar  to  the  shear  plates.  A  large  pressure 
gradient  builds  on  the  explosion  side  of  the  lower  face  and  a  shadow,  or  wake  forms 
behind  the  wedge  in  the  hollow  interior.  As  the  flow  progresses  past  the  end  of  the 
wedge  diffraction  begins  and  the  shadow  disappears.  The  main  difference  between  the 
shear  plates  and  the  wedge  obstacle  is  that,  due  to  its  size  and  geometry,  the  wedge  does  a 
better  job  of  lowering  the  pressure  in  its  wake.  The  wake,  however,  disappears  rapidly, 
due  to  diffraction,  after  the  shock  passes  the  end  of  the  obstacle.  The  orientation  of  the 
obstacle  is  not  as  critical  as  with  the  shear  plates.  The  figure  below  illustrates  these 
points 
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Figure  9  Flow  field  results  for  wedge  obstacle 

The  previous  figure  illustrates  the  effects  of  a  wedge  obstacle  on  the  flow.  The  geometry 
of  the  wedge  is  such  that  one  face  is  roughly  perpendicular  to  the  shock  and  a  large 
pressure  gradient  builds  on  this  face.  This  is  shown  by  the  closely  spaced  contour  lines  in 
the  lower  left  corner  between  grid  (4,4)  and  (7,2).  The  effects  of  this  pressure  gradient 
are  that,  the  wedge  absorbs  a  large  amount  of  energy,  the  flow  is  split  and  a  shadow  or 
wake  forms  inside  the  faces  of  the  wedge.  The  pressure  gradient  downstream  of  the 
wedge  is  'softened',  but  its  magnitude  is,  once  again,  similar  to  the  unimpeded  shock. 
Once  again,  following  passage  of  the  shock  past  the  end  of  the  wedge  diffraction  begins 
and  the  pressure  begins  to  fill  the  shadow.  This  result  is  not  unexpected. 


Since  the  orientation  of  the  wedge  is  not  as  critical  to  its  performance  as  that  of  the  shear 
plates,  wedges  would  be  more  a  practical  protection  than  the  other  obstacles  investigated. 
However,  due  to  the  low  viscosity  of  air  the  wedge  would,  in  all  likelihood  have  to  be 
incorporated  as  an  element  on  the  exterior  of  the  structure,  and  not  a  separate  element. 
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This  incorporation  would  lead  to  a  structure  with  wedge  shaped  faces  which,  it  is  already 
known  to  demonstrate  good  resistance  to  surface  loading  from  many  sources  including 
blasts,  waves,  running  water  etc.  Finally  this  form  of  protection  is  neither  inexpensive 
nor  aesthetically  pleasing. 


CHAPTER  6:  COST  ANALYSIS 


The  cost  analyses  for  the  obstacles  investigated  are  moot.  Since  the  pole  obstacle  does 
not  work  as  hoped,  unless  many  poles  are  placed  closely  together,  in  which  case  a  wall 
would  most  likely  be  cheaper,  I  will  present  a  synopsis  of  the  available  cost  analyses  for 
the  conventional  approach  to  blast  protection.  The  following  text  is  a  compilation  of 
others  work  and  is  not  to  be  considered  my  own.  These  are  my  words  but  not  my  ideas  or 
effort. 

Nuclear  Disasters  and  the  Built  Environment  contains  an  overview  of  the  results  of  a 
testing  completed  in  the  1950's  with  above  ground  nuclear  weapons.  The  results  show 
that  an  overpressure  of  5  psi  (34.4  kPa)  will  completely  demolish  a  conventional  wood 
frame  house  and  an  overpressure  of  1.7  psi  (11.7  kPa)  will  result  in  serious  damage.  A 
later  test  series  included  strengthened  wood  frame  houses.  The  first  house,  which  was 
subjected  to  4  psi  (27.6  kPa)  overpressure,  sustained  serious  structural  damage  but 
remained  standing.  The  second  house  was  exposed  to  2.6  psi  (17.9  kPa)  and  suffered 
damage  similar  to  the  unreinforced  house  exposed  to  1 .7  psi.  Also  included  in  the  second 
test  were  two  brick  houses.  These  were  found  to  be  of  the  same  strength  as  the  wood 
frame  houses.  One  brick  house  was  exposed  to  5  psi  (34.4  kPa)  and  sustained  damage 
similar  to  the  wood  frame  house.  The  second  house  received  an  overpressure  of  1.7  psi 
(11.7  kPa)  and  sustained  much  less  damage  than  the  wood  frame  house. 

The  strengthening  of  the  houses  for  the  test  resulted  in  a  5%  construction  premium  over 
unreinforced  construction. 

The  table  below  compiles  the  strengths  of  some  common  building  materials.  As  can  be 
seen  the  overpressure  required  for  failure  is  quite  low. 
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Table  1 1  Strength  of  common  building  materials 


Structural  element 

Overpressure  to  cause  failure 
(psi) 

Unreinforced  glass 

0.5-1.0 

Corrugated  steel  or  aluminum  siding 

1.0-2.0 

Brick  wall,  unreinforced 

3.0-10.0 

Standard  house  construction 

1.0-2.0 

Concrete  wall  panels  8"- 12"  thick 

1.5-5.5 

Protecting  Buildings  from  Bomb  Damage  contains  an  extensive  and  thorough  cost 
analysis  of  hardening  a  new  building  to  resist  exterior  explosions.  I  will  summarize  the 
results  of  the  analysis.  The  model  was  based  upon  a  5%  construction  premium  for 
hardening  the  building  and  also  assumed  a  minimum  10%  return  on  investment.  The 
construction  premium  is  the  additional  cost  that  would  be  incurred  if  the  blast  hardening 
were  to  be  included  in  the  construction  of  a  building.  The  construction  premium 
assumption  is  in  keeping  with  the  result  from  the  tests  completed  in  the  1950's.  The 
output  from  the  model  was  the  lease  premium  based  upon  these  assumptions.  These 
results  are  shown  in  the  table  below. 
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Table  12  Construction-lease  premium 


Assume  construction  premium 

Resulting  lease  premium 

3% 

2.69% 

5% 

3.46% 

7% 

4.23% 

The  lease  premium  in  the  table  above  is  the  additional  amount  that  would  have  to  be 
charged  lessees  in  a  blast  hardened  building.  This  lease  premium  would  recover  the 
construction  premium  at  the  1 0%  return  on  investment.  The  table  clearly  shows  that  the 
resulting  lease  premium  is  less  than  the  construction  premium.  Therefore  it  can  be 
economically  feasible  to  provide  blast  protection  in  new  construction,  if  the  developer  or 
customer  deems  it  necessary. 
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APPENDIX  A:  MACCORMACK  1-D  CODE  LISTING 


% 

%                                            Some  constants  needed  for  calculations  (UNITS) 

% 

%Courant  number 

CFL=input('Enter  the  Courant  number(0<CFL<  1  0.5  is  a  good  starting  point):  '); 

ifCFL>l. 

error('The  Courant  number  must  be  less  than  1 !') 

elseifCFL<=0. 

error('The  Courant  number  must  be  larger  than  zero!') 

end 
% 

%gamma  for  air 
g=1.4; 

gi=g-i; 

% 

%rho  =  initial  density 

(kg/mA3) 
rho=1.614; 
% 

%Pin  =  initial  pressure  (Pa) 

Pin=101325; 

% 

%T=initial  temperature 

(K) 
Tl=298; 
% 

%calculate  the  speed  of  sound  at  the  initial  pressure 

(m/s) 
c=(1.4*Pin/rho)A(l/2); 
% 

%h=  distance  over  which  to  perform  calculations  (m) 

h=input  ('Enter  the  distance  in  meters  over  which  to  perform  calculations:  '); 

ifh<=0. 

eiTor('You  must  enter  a  nonzero  distance!') 

end 
% 

%stop=  number  of  position  steps  per  distance 

stop=input  ('Enter  the  total  number  of  x  grid  points(  101  is  a  good  starting  point):  '); 
ifstop<=0. 
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error('You  must  enter  a  nonzero  number  of  x  grid  points!') 
end; 
% 

%calculate  the  x  step  size 

dx=h/(stop-l); 

% 

%                                                                   Print  output  of  values  on  screen 
% 

fprintfC W) 

fprintfC  Calculation  values 

\n') 

fprintf(The  Courant  number  is %1.2f\n',CFL) 

fprintf(The  distance  to  calculate  over  is %3.2f  m\n',h) 

fprintf(The  number  of  x  grid  points  is %4.0fvn',stop-l) 

fprintf('The  x  step  size  is %1.4f  m\n',dx) 

fprintfC W) 

forintfr  ***************************************************************** 

*\n') 

fprintfC  Conditions  ahead  of  the  shock 

V) 

fprintfCGamma  is %  1 . 1  f\n',g) 

fprintf(The  initial  pressure  is %3.3g  Pa\n',Pin) 

fprintf('The  initial  density  is %1.4f  kg/mA3\n',rho) 

fprintf(The  initial  temperature  is %3.3g  K\n',Tl) 

fprintf(The  speed  of  sound  is %3.4f  m/s\n',c) 

forintff * **************************************************************** 

*\n') 

% 

endt=input('Enter  the  number  of  time  steps  to  calculate  over:  "); 

ifendt<=0 

error('You  must  enter  a  number  of  time  steps  larger  than  zero!1) 

end 
% 

%define  vector  size 

oldp=[l:stop]*0.0; 

newp=[l:stop]*0.0; 

oldr=[l:stop]*0.0; 

newr=[l  :stop]*0.0; 

barr=[l:stop]*0.0; 

oldru=[l:stop]*0.0; 

newru=[l:stop]*0.0; 
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barru=[l:stop]*0.0; 

newu=[l:stop]*0.0; 

oldE=[l:stop]*0.0; 

newE=[l:stop]*0.0; 

barE=[l:stop]*0.0; 

loc=[l:stop]*0.0; 

% 

%  Convert  x  grid  points  to  actual  distances 

for  x=  1 :  stop 

loc(x)=(x-l)*dx; 
end 
% 

%  Calculate  Peak  staic  overpressure  Ps  based  on  charge  size 

%  and  range  to  charge  the  driving  force  for  the 

%  calculations  to  follow 

% 

%  r=range  to  charge 

(m) 
r=input('Enter  the  range  to  the  charge  in  meters:  '); 

ifr<=0 

error('You  must  enter  a  range  larger  than  zero!') 

end 
% 

%                     Print  out  table  of  explosive  types 
%                                            T=type  of  explosive 
fprintfC W) 

fprintfC  Choose  the  explosive  using  the  table  below\t\n') 

fprintf('Compound  B l\n') 

fprintf('RDX  (Cyclonite) 2V) 

fprintf('HMX 3\n') 

fprintf('Nitroglycerin  (liquid) 4\n') 

fprintf(TNT 5V) 

fprintf('Blasting  Gelatin 6\n') 

fprintf('60  percent  Nitroglycerin  dynamite 7\n') 

fprintf('Semtex 8\n') 

fprintfC \n') 

% 

%                      Ask  for  whiat  type  of  explosive  to  use 
% 

T=input('enter  the  type  of  explosive  to  use:  '); 
if  T==l 

S=1.148; 


51 

elseif  T==2 

S=1.185; 
elseif  T==3 

S=1.256; 
elseif  T==4 

S=1.481; 
elseif  T==5 

S=1.000; 
elseif  T==6 

S=1.000; 
elseif  T==7 

S=0.600; 
elseif  T==8 

S=1.250; 
else 

eiTor('You  must  enter  an  explosive  type!') 
end 

fprintf(The  TNT  conversion  factor  for  this  explosive  is:  %1.3f\n',S) 
% 

%                     Ask  what  charge  mass  to  use  (kg  of  TNT) 

% 

M=input('Enter  the  mass  of  the  charge  in  kilograms:  '); 
ifM<0 

error('You  must  enter  a  mass  larger  than  zero!') 
end 

W=M*S; 

% 

%  calculate  peak  static  overpressure  (Bar) 

z=r/(WA(l/3)); 
Ps=(6.7/(zA3)); 

if  Ps<10 

Ps=((0.975/z)+(  1 .455/(zA2))+(5. 85/(zA3)))-0.0 1 9; 

else  Ps=(6.7/(zA3)); 

end 

if  Ps<0.1 

Ps=Ps; 

end 

% 

%  Convert  Ps(the  static  overpressure)  from  Bar  (Pa) 

%  to  Pascals  from  overpressure  to  absolute  pressure 
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%                      and  apply  conversion  factor  for  ground  burst 
% 

Ps=1.8*((Ps*100000)+Pin); 


% 

%                                                                   Print  output  of  values  on  screen 
fprintfC V) 

fprintf('  Explosive  values 

\n') 

fprintf(The  range  to  the  charge  from  the  x=0  point  is %4.2f  m\n',r) 

fprintf(The  mass  of  the  charge  in  TNT  equivalents  is %4.3f  kg\n',W) 

fprintf('The  static  overpressure  of  the  charge  is %1.3g  Pa\n',Ps-Pin) 

fprintfC V) 

% 

%                                                                      Calculate  intial  velocity  using  Ps 
% 

Mach=sqrt(((Ps/Pin)- 1  )*((g+ 1  )/(2*g))+ 1 ); 

Vin=c*(2/(g+l)*(Mach-(l/Mach))); 

%  Print  output  of  values  on  screen 

fprintf('The  initial  wave  velocity  is %4.4f  m/s\n',Vin) 

fprintf('The  shock  wave  Mach  number  is %2.4f\n',Mach) 

fprintfC V) 

% 

%                                                                   Calculate  values  behind  the  shock 
% 

rho2=rho/(  1  -(2/(g+ 1  ))*(  1  -( 1  /MachA2))); 

T2=Ps/(rho2*(83 14.51/29)); 

c2=(cA2*T2/Tl)A0.5; 

forintff '  *  **************************************************************** 

*W) 

fprintfC  Conditions  behind  the  shock 

W) 

fprintf(The  pressure  is %1.3g  Pa\n',Ps) 

fprintf(The  density  is %1.4f  kg/mA3\n',rho2) 

fprintf(The  temperature  is %3.3g  K\n',T2) 

fprintf(The  speed  of  sound  is %3.4f  m/s\n',c2) 

forintff  ****************************************************************** 

*\n') 

% 

%                                                                   Define  initial  conditions 
% 


53 


for  x=l:stop 

oldp(x)=Pin; 
oldr(x)=rho; 
oldru(x)=0; 
oldE(x)=Pin/gl; 

end 
% 

%                     Define  boundary  conditions  at  grid  point  x=l 
% 

oldp(l)=Ps; 

oldr(  1  )=rho/(  1  -<2/(g+ 1  ))*(  1  -( 1  /MachA2))); 

oldru(l)=oldr(l)*Vin; 

oldE(l)=(Ps/gl)+(0.5*oldru(l)A2/oldr(l)); 


% 

%  Calculate  the  new  speed  of  sound  based  on  new  pressure 

%  and  density  values 

% 

newc=sqrt(  1 .4*oldp(x)/oldr(x)); 
% 

%                                                                   Calculate  the  initial  time  step 
% 

dt=CFL*dx/(newc+oldru(  1  )/oldr(  1)); 

fprintf(The  inital  time  step  is %1.5es\n',dt) 

% 

%                                                                   Initialize  total  time  counter 
% 

tt=0; 

% 

%  Define  predictor  values  at  x  equal  1 

% 

barr(  1  )=oldr(  1  )-(dt)*(oldru(2)-oldru(  1 )); 
barru(  1  )=oldru(  1  )-(dt)*(((oldru(2)A2/oldr(2))-. . . 

(oldru(  1  )A2/oldr(  1  )))+(oldp(2)-oldp(  1 ))); 
barE(  1  )=oldE(  1  )-(dt)*((oldE(2)+oldp(2))*oldru(2)/.. . 

oldr(2)-(oldE(  1  )+oldp(  1  ))*oldru(  1  )/oldr(  1 ) ); 

% 

%  Adjust  x  grid  end  point  for  finite  differencing 

stop=stop- 1 ; 
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%                      Finite  difference  for  density,  momentum  and  energy 
% 


for  t=  1 :  endt 

dtmin=l.; 

CFLmin=CFL; 

dtdx=dt/(dx); 
%Calculate  total  time 

tt=tt+dt; 

for  x=2:  stop 

% 

%  PREDICTOR 

% 

%  Solve  as  a  function  of  time  and  position  for  increasing 

%  time 

% 

%  Density(r) 

barr(x)=oldr(x)-(dtdx)*(oldru(x+ 1  )-oldm(x)); 
% 

%  Momentum(ru) 

barru(x)=oldru(x)-(dtdx)*((oldru(x+ 1 )  A2/oldr(x+ 1 ))-. . . 

(oldru(x)A2/oldr(x))+(oldp(x+ 1  )-oldp(x))); 
% 

%  Energy(E) 

barE(x)=oldE(x)-(dtdx)*((oldE(x+ 1  )+oldp(x+ 1  ))*. . . 
(oldru(x+l)/oldr(x+l))-(oldE(x)+oldp(x))*... 
(oldru(x)/oldr(x))); 
% 

%  CORRECTOR 

% 

%  Density(r) 

newr(x)=(barr(x)+oldr(x))/2-(dtdx/2)*(barru(x)-barru(x- 1 )); 

%  Left  wall  boundary  conditions 

newr(stop+ 1  )=newr(stop- 1 ); 
newr(stop)=newr(stop-l); 
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% 

%  Momentum(ru) 

newru(x)=(barru(x)+oldru(x))/2-(dtdx/2)*((banu(x)A2/... 
barr(x))-(barru(x- 1 )  A2/barr(x- 1  ))+(g  1  *(barE(x)-. . . 
(.5*(ban-u(x)A2/barr(x)))))-(gl*(barE(x-l)-(5*... 
(barru(x- 1  )A2/baiT(x- 1 )))))); 

%  Left  wall  boundary  conditions 

newru(stop)=0; 
newru(stop+ 1  )=0; 

% 

%  Energy(e) 

newE(x)=(barE(x)+oldE(x))/2-(dtdx/2)*((barru(x)/... 
baiT(x))*(barE(x)+(gl*(barE(x)-(.5*(ban-u(x)A2/... 
baiT(x))))))-((barru(x- 1  )/barr(x- 1  ))*(barE(x- 1 )+. . . 
(gl*(barE(x-l)-(.5*(barru(x-l)A2/barr(x-l )))))))); 

% 

%  Calculate  pressure  using  the  continuity  equation  for  the 

%  ideal  gas  law  and  the  new  values  for  energy,  momentum 

%  and  density 

newp(x)=gl*(newE(x)-(.5*(newru(x)A2/newr(x)))); 

newp(stop+l)=newp(stop-l); 


% 

%  Calculate  velocity  using  the  new  values  for  momentum 

%  and  density 

newu(x)=newru(x)/newr(x); 
% 

%  Check  stability 

%  Calculate  the  new  speed  of  sound 

newc=sqrt(g*newp(x)/newr(x)); 
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% 

%  Calculate  the  new  time  step  dt 

dtnew=dx/(newu(x)+newc); 
if  dtnew  <  dtmin 

dtmin=dtnew; 
end 

CFLnew=(dtmin/dx)*((newc+oldru(x)/oldr(x))); 
if  CFLnew<CFLmin 

CFLmin=CFLnew; 
end 
% 

%  Return  for  next  x 

end 

% 

%  Apply  numerical  boundary  scheme  to  calculate  all  values 

%  atx=l  andx=stop+l 

newr(l)=oldr(l); 

newE(l)=oldE(l); 

newru(l)=oldru(l); 

newu(l)=oldru(l)/oldr(l); 

newp(l)=oldp(l); 

barr(  1  )=oldr(  1  )-(dt)*(oldru(2)-oldru(  1 )); 

bai-ru(  1  )=oldru(  1  )-(dt)*((oldru(2)A2/oldr(2)-. . . 

oldru(l)A2/oldr(l))+(oldp(2)-oldp(2))); 
barE(  1  )=oldE(  1  )-(dt)*((oldE(2)+oldp(2))*oldru(2)/.. . 

oldr(2)  -(oldE(l)+oldp(l))*oldru(l)/oldr(l) ), 
newE(stop+ 1  )=oldE(stop+ 1 ); 
newu(stop+ 1  )=newru(stop+ 1  )/newr(stop+ 1 ); 


% 

%  Update  time  step  and  CFL(Courant  number) 

dt=dtmin; 

CFL=CFLmin; 
fprintf('CFL=  %1.3g\tdeltat=  %1.5e  s\tTotal  time=  %1.5e  s\n',CFL,dt,tt) 


%- 
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%  Output  the  values  for  pressure,  density  and  velocity 

%  using  graphs 

% 

%  Pressure 

subplot(3,l,l),  plot(loc,newp),  ylabel('Absolute  Pressure  (Pa)') 

% 

%  Density 

subplot(3,l,2),  plot(loc,newr),  ylabel('Density  (kg/mA3)') 

% 

%  Velocity 

subplot(3,l,3),  plot(loc,newu),  ylabel('Velocity  (m/s)'),... 
xlabel('Distance  (m)') 

% 

pause(.OOOl) 

% 

%Replace  old  values  with  new  values  and  begin  next  time  step 

oldp=newp; 

oldr=newr; 

oldru=newru; 

oldE=newE, 

% 

%  Return  for  next  t 

end 

% 


APPENDIX  B:  MACCORMACK  QUASI  2-D  CODE  LISTING 
clear  all 
% 

%  Some  constants  needed  for  calculations  (UNITS) 

%Courant  number 

CFL=input('Enter  the  Courant  number(0<CFL<l  0.5  is  a  good  starting  point):  '); 

ifCFL>l. 

error(The  Courant  number  must  be  less  than  1 !') 

elseifCFL<=0. 

eiTor(The  courant  number  must  be  larger  than  zero!') 

end 
% 

%gamma  for  air 
g=14; 

gi=g-i; 

% 

%rho  =  initial  density 

(kg/mA3) 
rho=1.614; 

% 

%Pin  =  initial  pressure  (Pa) 

Pin=101325; 

% 

%T=initial  temperature 

(K) 
Tl=298; 

% 

%calculate  the  speed  of  sound  at  the  initial  pressure  (m/s) 

c=(1.4*Pin/rho)A(l/2); 

% 

%h=  distance  over  which  to  perform  calculations  (m) 

h=input  ('Enter  the  distance  in  meters  over  which  to  perform  calculations:  '); 

ifh<=0. 

error('You  must  enter  a  nonzero  distance!') 

end 
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% 

%stop=  number  of  position  steps  per  distance 

stop=input  ('Enter  the  total  number  of  grid  points  in  the  x,y  directions(3 1  is  a  good 

starting  point):  '); 

ifstop<=0. 

error('You  must  enter  a  nonzero  number  of  x  grid  points!') 

end; 

% 


%calculate  the  x  and  y  step  size 

dx=h/(stop-l); 

dy=h/(stop-l); 


%- 


% 

fprintfO- 

fprintf(' 


Print  output  of  values  on  screen 
\n') 

Calculation  values 


\n') 


fprintf('The  Courant  number  is %1.2f\n',CFL) 

fprintf('The  x  distance  to  calculate  over  is %3.2f  m\n',h) 

fprintf('The  y  distance  to  calculate  over  is %3.2f  m\n',h) 

fprintf('The  number  of  x  grid  points  is %4.0f\n',stop-l) 

fprintf(The  number  of  y  grid  points  is %4.0f\n',stop-l) 

fprintf(The  x  step  size  is %1.4f  m\n',dx) 

fprmtf(The  y  step  size  is %1.4f  m\n',dx) 

fprintf(' \n') 

*\n*) 

fprintf('  Conditions  ahead  of  the  shock 

\n') 

fprintf('Gamma  is %l.lf\n',g) 

fprintf('The  initial  pressure  is %3.3g  Pa\n',Pin) 

fprintf('The  initial  density  is %1.4f  kg/mA3\n',rho) 

fprintf('The  initial  temperature  is %3.3g  K\n',Tl) 

fprintf(The  speed  of  sound  is %3.4f  m/s\n',c) 

*\n') 


%- 


endt=input('Enter  the  number  of  time  steps  to  calculate  over:  '); 
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ifendt<=0 

error('You  must  enter  a  number  of  time  steps  larger  than  zero!') 

end 

% 

%  define  vector  size 

for  x=l:stop 

for  y=l:stop 

op(x,y)=0.0; 

np(x,y)=0.0; 

Olr(x,y)=0.0; 

nr(x,y)=0.0; 

br(x,y)=0.0; 

Olru(x,y)=0.0; 

nru(x,y)=0.0; 

Olrv(x,y)=0.0; 

nrv(x,y)=0.0; 

bru(x,y)=0.0; 

nu(x,y)=0.0; 

brv(x,y)=0.0; 

nv(x,y)=0.0; 

oE(x,y)=0.0; 

nE(x,y)=0.0; 

bE(x,y)=0.0; 


end 
end 


% 

%                     Calculate  Peak  static  overpressure  Ps  based  on  charge 
%                       size  and  range  to  charge  the  driving  force  for  the 
%                       calculations  to  follow 
% 

%  r=range  to  charge 

(m) 
% 

r=input(' Enter  the  range  to  the  charge  in  meters:  '); 
ifr<=0 

error('You  must  enter  a  range  larger  than  zero!') 
end 
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%  Print  out  table  of  explosive  types 

%  T=type  of  explosive 

fprintfC \n') 

fprintf('  Choose  the  explosive  using  the  table  below\t\n') 

fprintf('Compound  B l\n') 

fprintf('RDX  (Cyclonite) 2\n') 

fprintf('HMX 3\n') 

fprintf('Nitroglycerin  (liquid) 4\n') 

fprintf('TNT 5\n*) 

fprintf('Blasting  Gelatin 6\n') 

fprintf('60  percent  Nitroglycerin  dynamite 7\n') 

fprintf('Semtex 8V) 

fprintfC \n') 

% 

%                      Ask  for  what  type  of  explosive  to  use 
% 

T=input('enter  the  type  of  explosive  to  use:  '); 
ifT=l 

S=1.148; 
elseif  T==2 

S=1.185; 
elseif  T==3 

S=1.256; 
elseif  T==4 

S=1.481; 
elseif  T==5 

S=1.000; 
elseif  T==6 

S=1.000; 
elseif  T==7 

S=0.600; 
elseif  T==8 

S=1.250; 
else 

error('You  must  enter  an  explosive  type!') 
end 
fprintf(The  TNT  conversion  factor  for  this  explosive  is:  %1 .3f\n',S) 

% 
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%  Ask  for  mass  of  charge 

(kg) 
% 

M=input('Enter  the  mass  of  the  charge  in  kilograms:  '); 

ifM<0 

error('You  must  enter  a  mass  larger  than  zero!') 

end 
W=M*S; 

% 

%                                            Calculate  peak  static  oveipressure  (B) 

% 

z=r/(WA(l/3)); 
Ps=(6.7/(zA3)); 

if  Ps<10 

Ps=((0.975/z)+(1.455/(zA2))+(5.85/(zA3)))-0.019; 

elsePs=(6.7/(zA3)); 

end 

if  Ps<0.1 

Ps=Ps; 

end 

% 

%  Convert  Ps(the  static  oveipressure)  from  Bar  (Pa) 

%  to  Pascals  from  overpressure  to  absolute 

%  pressure  and  apply  conversion  factor  for 

%  ground  burst 

% 

Ps=1.8*((Ps*100000)+Pin); 


% 

%                                                                      Print  output  of  values  on  screen 
% 

fprintfC \n') 

fprintf('  Explosive  values 

V) 
fprintf(The  range  to  the  charge  from  the  x=0  point  is %4.2f  m\n',r) 
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fprintf(The  mass  of  the  charge  in  TNT  equivalents  is %4.3f  kg\n',W) 

fprintf(The  static  overpressure  of  the  charge  is %1.3g  Pa\n',Ps-Pin) 

fprintfC \n') 

% 

%                                                                      Calculate  intial  velocity  using  Ps 
% 

Mach=sqrt(((Ps/Pin)- 1  )*((g+ 1  )/(2*g))+ 1 ); 

Vin=c*(2/(g+l)*(Mach-(l/Mach))); 

%  Print  output  of  values  on  screen 

fprintf(The  initial  wave  velocity  is %4.4f  m/s\n',Vin) 

fprintf(The  shock  wave  Mach  number  is %2.4f\n',Mach) 

fprintfC V) 

% 

%                                                                   Calculate  values  behind  the  shock 
% 

rho2=rho/(  1  -(2/(g+ 1  ))*(  1  -( 1  /MachA2))); 

T2=Ps/(rho2  *(83 14.51  /29)); 

c2=(cA2*T2/Tl)A0.5; 

forintfr  *  **************************************************************** 

*W) 

fprintfC  Conditions  behind  the  shock 

W) 

fprintf(The  pressure  is %1.3g  Pa\n',Ps) 

fprintf(The  density  is %1.4f  kg/mA3V,rho2) 

fprintf(The  temperature  is %3.3g  K\n*,T2) 

fprintf(The  speed  of  sound  is %3.4f  m/s\n',c2) 

forintff  '***************************************************************** 

*\n') 

% 

% 

%                                                                      Define  initial  conditions 
% 

for  x=l:stop 

for  y=l:stop 

op(x,y)=Pin; 

01r(x,y)=rho; 

Olru(x,y)=0; 
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Olrv(x,y)=0; 
oE(x,y)=Pin/gl; 
end 

end 

% 

%                      Define  boundary  conditions  at  grid  point  x=l  y=l 
% 

y=i; 

for  x=l:stop 

op(x,y)=Ps; 

01r(x,y)=rho/(  1  -(2/(g+ 1 ))  *(  1  -( 1  /MachA2))); 
Olru(x,y)=0; 
01rv(x,y)=01r(x,y)*Vin; 

oE(x,y)=(Ps/g  1  )+(0. 5  *(01ru(x,y)A2+01rv(x,y)A2)/. 
01r(x,y)); 

end 


% 

%  Calculate  the  new  speed  of  sound  based  on  new  pressure 

%  and  density  values 

% 

nc=sqrt(g*op(x,y)/01r(x,y)); 
% 

%                                                                   Calculate  the  initial  time  step 
% 

dt=CFL*(  l/(abs(01ru(  1 , 1  )/01r(  1 , 1  ))/dx+abs(01rv(  1 , 1  )/01r(  1 , 1 ))/. . . 
dy+nc*(  1  /dxA2+ 1  /dyA2)  A(0. 5))); 


fprintf(The  inital  time  step  is %1.5e  s\n',dt) 

% 

%                                                                   Initialize  total  time  counter 
% 

tt=0; 
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% 

%                      Finite  difference  for  density,  momentum  and  energy 
% 

% 

%                     Adjust  x  and  y  grid  end  point  for  finite  differencing 
% 

stop=stop- 1 ; 


for  t=  1 :  endt 

% 

%  NBS  for  column  1  and  row  1 

% 

dtmin=l.; 

dtdx=dt/dx; 
dtdy=dt/dy; 

% 

%                                                                      Calculate  total  time 
% 

tt=tt+dt; 

% 

for  x=l:  stop 
for  y=l:  stop 

% 

%  PREDICTOR 

% 

%  Solve  as  a  function  of  time  and  position  for 

%  increasing  time 
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%- 

% 

%- 


Density(r) 


yo- 
yo 
%- 


br(x,y)=01r(x,y)-... 
(dt*(((01ru(x+ 1  ,y)-01ru(x,y))/dx)+. . . 

((01rv(x,y+l)-01rv(x,y))/dy))); 

br(x,stop+ 1  )=br(x,stop); 
br(stop+l  ,y)=br(stop,y); 
br(stop+ 1  ,stop+ 1  )=br(stop,stop); 


Momentum(ru) 


%- 

% 

%- 


bru(x,y)=01ru(x,y)-dt*. . . 
((1/dx*... 

(01ru(x+ 1  ,y)A2/01r(x+ 1  ,y )-. . . 
01ru(x,y)A2  /01r(x,y)  +... 
(op(x+l,y)-op(x,y)        )))+... 
(1/dy*... 

(01rv(x,y+ 1  )*01ru(x,y+ 1  )/01r(x,y+ 1  )- 
01rv(x,y)  *01ru(x,y)  /01r(x,y)    ))); 

bru(x,stop+ 1  )=bru(x,stop); 
bru(stop+ 1  ,y)=bru(stop,y); 
bru(stop+ 1  ,stop+ 1  )=bru(stop,stop); 


Momentum(rv) 


brv(x,y)=01rv(x,y)-dt*. . . 
((1/dx*... 

(01ru(x+ 1  ,y)*01rv(x+ 1  ,y)/01r(x+ 1  ,y)- 
01ru(x,y)  *01rv(x,y)  /01r(x,y)  ))+... 
(1/dy*... 

(01rv(x,y+l)A2/01r(x,y+l)-... 
01rv(x,y)A2  /01r(x,y)+... 
(op(x,y+l)-op(x,y)       )))); 

brv(stop+l  ,y)=brv(stop,y); 
brv(x,stop+ 1  )=brv(x,stop); 
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brv(stop+ 1  ,stop+ 1  )=brv(stop,stop); 
% 

%  Energy  (e) 

% 

bE(x,y)=oE(x,y)-dt*... 
((1/dx*... 

((oE(x+l,y)+op(x+l,y))*01m(x+l,y)/01r(x+l,y)  -... 
(oE(x,y)+op(x,y)     )*01ru(x,y)  /01r(x,y)    )+... 
1/dy*... 
((oE(x,y+l)+op(x,y+l))*01rv(x,y+l)/01r(x,y+l)  -... 
(oE(x,y)+0p(x,y)     )*01rv(x,y)  /01r(x,y)    ))); 

bE(stop+ 1  ,y)=bE(stop,y); 

bE(x,stop+ 1  )=bE(x,stop); 

bE(stop+ 1  ,stop+ 1  )=bE(stop,stop); 
% 

%  return  for  next  x 

end 

%  return  for  next  y 

end 

% 

nr(:,l)=01r(:,l); 
nru(:,l)=01ru(:,l); 
nrv(:,l)=01rv(:,l); 
nE(:,l)=oE(:,l); 

for  x=2:stop+l 

for  y=2:stop+l 

% 

%  CORRECTOR 

% 

%  Density(r) 

% 

nr(x,yHbr(x,y)+01r(x,y))/2-... 
(  dtdx/2*(bru(x,y)-bru(x-l,y))+... 
dtdy/2*(brv(x,y)-brv(x,y-l))); 
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nr(l,y)=nr(2,y); 

%  Left  wall  boundary  conditions 

% 

%  Momentum(ru) 

% 

Px=gl*(bE(x,y)-0.5*((bru(x,y)A2/br(x,y))+(bi-v(x,y)A2/br(x,y)))); 

Pxm=g  1  *(bE(x- 1  ,y)-0. 5*((bru(x- 1  ,y)A2/br(x- 1  ,y))+(brv(x- 1  ,y)A2/br(x- 1  ,y)))); 

nru(x,y)=(bru(x,y)+01m(x,y))/2-... 
((dtdx/2*... 

(( bru(x,y)  *bru(x,y)  /br(x,y)  +Px)-... 
(  bru(x- 1  ,y)  *bru(x- 1  ,y)/br(x- 1  ,y)+Pxm)))+. . . 
(dtdy/2*... 

(  bru(x,y)  *brv(x,y)  /br(x,y)-... 
bru(x,y-l)*brv(x,y-l)/br(x,y-l)      ))) ; 

% 


nru(l,y)=nru(2,y); 

%  Momentum(rv) 

% 

Py=gl*(bE(x,y)-0.5*((bru(x,y)A2/br(x,y))+(brv(x,y)A2/br(x,y)))); 

Pym=g  1  *(bE(x,y- 1  )-0. 5  *((bru(x,y- 1  )A2/br(x,y- 1  ))+(brv(x,y- 1  )A2/br(x,y- 1 )))); 

nrv(x,y)=(brv(x,y)+01rv(x,y))/2-... 
((dtdx/2*... 
(  bru(x,y)  *bi^v(x,y)  /br(x,y) 

bru(x- 1  ,y)*brv(x- 1  ,y)/br(x- 1  ,y) ))+... 
(dtdy/2*... 
( (brv(x,y)*brv(x,y)/br(x,y)+      Py  )-... 

(brv(x,y- 1 )  *brv(x,y- 1  )/br(x,y- 1  )+Pym)))); 

% 


nrv(l,y)=nrv(2,y); 
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%  Left  wall  boundary  conditions 

% 

%  Energy  (E) 

% 

nE(x,y)=(bE(x,y)+oE(x,y))/2-... 
((dtdx/2*... 

((bru(x,y)  /br(x,y)  *(bE(x,y)  +Px  ))    -... 
(bru(x- 1  ,y)/br(x-l  ,y)*(bE(x- 1  ,y)+Pxm))))  +. . . 
(dtdy/2*... 

((bi-v(x,y)  /br(x,y)  *(bE(x,y)+Py    ))    -... 
(brv(x,y- 1  )/br(x,y- 1  )*(bE(x,y- 1  )+Pym))))    ); 
% 

%— 

%  Right  wall  boundary  conditions 

% 

nr(x,stop+l)=nr(x,stop-l); 
%  nr(x,stop)=nr(x,stop-l); 

nrv(x,stop)=0; 
nrv(x,stop+l)=0; 
nE(x,stop+ 1  )=oE(x,stop+ 1 ); 

nE(l,y)=nE(2,y); 
% 

%return  for  next  x 
end 
% 

%return  for  next  y 
end 
% 

for  y=l:stop+l 
for  x=l:stop+l 
% 

%                      Calculate  pressure  using  the  continuity  equation  for 
%                       the  ideal  gas  law  and  the  new  values  for  energy, 
%                       momentum  and  density 
% 
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np(x,y)=(gl)*(nE(x,y)-(.5*((nm(x,y)A2/. 
nr(x,y))+(nrv(x,y)A2/nr(x,y))))); 

np(x,stop+ 1  )=np(x,stop- 1 ), 


% 

%  Calculate  velocity  using  the  new  values  for  momentum 

%  and  density 

% 

nu(x,y)=nru(x,y)/nr(x,y); 
nv(x,y)=nrv(x,y)/nr(x,y); 


% 

%  Check  stability 

%                                                                   Calculate  the  new  speed  of  sound 
% 

nc=sqrt(g*np(x,y)/nr(x,y)); 
% 

% 

%                                                                      Calculate  the  new  time  step  dt 
% 


dtn=CFL*(  1  /(abs(nu(x,y  ))/dx+abs(nv(x,y))/dy+. 

nc*(  l/dxA2+l/dyA2)A(0.5))); 
if  dtn  <  dtmin 

dtmin=dtn; 
%  end  dtmin  if  statement 

end 
% 

% 

% 

%return  for  next  x 

end 
% 

%return  for  next  y 
end 
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%- 


% 

%  Update  time  step  and  CFL(Courant  number) 

% 

dt=dtmin; 
fprintf('delta  t=  %1.5e  s\tTotal  time=  %1.5e  s\n',dt,tt) 


% 

%  Output  the  values  for  pressure,  density  and  velocity 

%  using  graphs 

% 

%  Pressure 


colormap(cool) 
surf(np) 
hold  on 

pcolor((np)-Pin) 
title('Absolute  Pressure') 
xlabel('X  grid  points') 
ylabel('Y  grid  points') 
zlabel('Pascals') 

pause(.Ol) 

hold  off 

% 

%  Replace  old  values  with  new  values  and  begin 

%  next  time  step 

% 

op=np; 

01r=nr; 

01ru=nru; 

01rv=nrv; 

oE=nE; 

% 
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%  Return  for  next  t 

end 

% 


APPENDIX  C:  MACCORMACK  2-D  CODE  LISTING 


clear  all 

o/0 


%  Some  constants  needed  for  calculations  (UNITS) 

%Courant  number 

CFL=input('Enter  the  Courant  number(0<CFL<l  0.5  is  a  good  starting  point):  '); 

ifCFL>l. 

eiTor('The  Courant  number  must  be  less  than  1 !') 

elseifCFL<=0. 

eiTor(The  courant  number  must  be  larger  than  zero!') 

end 
% 

%gamma  for  air 
g=1.4; 

gi=g-i; 

% 

%rho  =  initial  density 

(kg/mA3) 
rho=1.614; 

% 

%Pin  =  initial  pressure  (Pa) 

Pin=101325; 

% 

%T=initial  temperature 

(K) 
Tl=298; 

% 

%calculate  the  speed  of  sound  at  the  initial  pressure  (m/s) 

c=(1.4*Pin/rho)A(l/2); 

% 

%h=  distance  over  which  to  perform  calculations  (m) 

h=input  ('Enter  the  distance  in  meters  over  which  to  perform  calculations:  '); 

ifh<=0. 

error('You  must  enter  a  nonzero  distance!') 

end 
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% 

%stop=  number  of  position  steps  per  distance 

stop=input  ('Enter  the  total  number  of  grid  points  in  the  x,y  directions(3 1  is  a  good 

starting  point):  '); 

ifstop<=0. 

error('You  must  enter  a  nonzero  number  of  x  grid  points!') 

end; 


% 

%calculate  the  x  and  y  step  size 

dx=h/(stop-l); 

dy=h/(stop-l); 


%- 


%  Print  output  of  values  on  screen 

fprintfC \n') 

fprintfC  Calculation  values 

\n') 

fprintf(The  Courant  number  is %1.2f\n',CFL) 

fprintf('The  x  distance  to  calculate  over  is %3.2f  m\n',h) 

fprintf('The  y  distance  to  calculate  over  is %3.2fm\n',h) 

fprintf('The  number  of  x  grid  points  is %4.0f\n',stop-l) 

fprintf('The  number  of  y  grid  points  is %4.0f\n',stop-l) 

fprintf('The  x  step  size  is %1.4f  m\n',dx) 

fprintf(The  y  step  size  is %1.4f  m\n',dx) 

fprintf(' \n') 

forintfr  *  **************************************************************** 

*\n') 

fprintf('  Conditions  ahead  of  the  shock 

\n') 

fprintf('Gamma  is %l.lf\n',g) 

fprintf('The  initial  pressure  is %3.3g  Pa\n',Pin) 

fprintf('The  initial  density  is %1.4f  kg/mA3\n',rho) 

fprintf(The  initial  temperature  is %3.3g  K\n',Tl) 

fprintf(The  speed  of  sound  is %3.4f  m/s\n',c) 

forintfr  ****************************************************************  * 

*\n') 

% 
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endt-input('Enter  the  number  of  time  steps  to  calculate  over:  '); 
ifendt<=0 

eiTor('You  must  enter  a  number  of  time  steps  larger  than  zero!') 
end 

% 

%  define  vector  size 

for  x=l:stop 

for  y=l:stop 

op(x,y)=0.0; 

np(x,y)=0.0; 

Olr(x,y)=0.0; 

nr(x,y)=0.0; 

br(x,y)=0.0; 

Olru(x,y)=0.0; 

nru(x,y)=0.0; 

Olrv(x,y)=0.0; 

nrv(x,y)=0.0; 

bru(x,y)=0.0; 

nu(x,y)=0.0; 

brv(x,y)=0.0; 

nv(x,y)=0.0; 

oE(x,y)=0.0; 

nE(x,y)=0.0; 

bE(x,y)=0.0; 


end 
end 


% 

%                     Calculate  Peak  static  overpressure  Ps  based  on  charge 
%                       size  and  range  to  charge  the  driving  force  for  the 
%                      calculations  to  follow 
% 

%  r=range  to  charge 

(m) 
% 

r=input('Enter  the  range  to  the  charge  in  meters:  '); 
ifr<=0 
error('You  must  enter  a  range  larger  than  zero!1) 
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end 

%                      Print  out  table  of  explosive  types 
%                                            T=type  of  explosive 
fprintfC \n') 

fprintf('  Choose  the  explosive  using  the  table  below\t\n') 

fprintf('Compound  B IV) 

fprintf('RDX  (Cyclonite) 2\n') 

fprintf('HMX 3\n') 

fprintf('Nitroglycerin  (liquid) 4\n') 

fprintf(TNT 5\n') 

fprintf('Blasting  Gelatin 6\n') 

fprintf('60  percent  Nitroglycerin  dynamite 7V) 

fprintf('Semtex 8\n*) 

fprintfC \n') 

% 

%                     Ask  for  what  type  of  explosive  to  use 
% 

T=input('enter  the  type  of  explosive  to  use:  '); 
ifT==l 

S=1.148; 
elseif  T==2 

S=1.185; 
elseif  T==3 

S=1.256; 
elseif  T==4 

S=1.481; 
elseif  T==5 

S=1.000; 
elseif  T==6 

S=1.000; 
elseif  T==7 

S=0.600; 
elseif  T==8 

S=1.250; 
else 

eiTor('You  must  enter  an  explosive  type!') 
end 
fprintf(The  TNT  conversion  factor  for  this  explosive  is:  %1.3f\n',S) 

% 
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%  Ask  for  mass  of  charge 

(kg) 
% 


M=input('Enter  the  mass  of  the  charge  in  kilograms:  '); 

ifM<0 

error('You  must  enter  a  mass  larger  than  zero!') 

end 
W=M*S; 

% 

%                                            Calculate  peak  static  overpressure  (B) 

% 

z=r/(WA(l/3)); 
Ps=(6.7/(zA3)); 

if  Ps<10 

Ps=((0.975/z)+(  1 .455/(zA2))+(5. 85/(zA3)))-0.0 1 9; 

else  Ps=(6.7/(zA3)); 

end 

if  Ps<0.1 

Ps=Ps; 

end 

% 

%  Convert  Ps(the  static  overpressure)  from  Bar  (Pa) 

%  to  Pascals  from  overpressure  to  absolute 

%  pressureand  apply  conversion  factor  for 

%  ground  burst 

% 

Ps=1.8*((Ps*100000)+Pin); 


% 

%                                                                      Print  output  of  values  on  screen 
% 

fprintfC \n') 

fprintf('  Explosive  values 

\n') 
fprintf(The  range  to  the  charge  from  the  x=0  point  is %4.2f  m\n',r) 
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fprintf(The  mass  of  the  charge  in  TNT  equivalents  is %4.3f  kg\n',W) 

fprintf('The  static  overpressure  of  the  charge  is %1.3g  Pa\n',Ps-Pin) 

fprintf(' W) 

% 

%                                                                      Calculate  intial  velocity  using  Ps 
% 

Mach=sqrt(((Ps/Pin)- 1  )*((g+ 1  )/(2*g))+ 1 ); 

Vin=c*(2/(g+ 1 )  *(Mach-(  1  /Mach))); 

%  Print  output  of  values  on  screen 

fprintf('The  initial  wave  velocity  is %4.4f  m/s\n',Vin) 

fprintf(The  Shockwave  Mach  number  is %2.4f\n',Mach) 

fprintf(' \n') 

% 

%                                                                      Calculate  values  behind  the  shock 
% 

rho2=rho/(  l-(2/(g+l))*(  1  -( l/MachA2))); 
T2=Ps/(rho2  *(83 14.51  /29)); 
c2=(cA2*T2/Tl)A0.5; 

*\n') 

fprintf('  Conditions  behind  the  shock 

W) 

fprintf(The  pressure  is %1.3g  Pa\n',Ps) 

fprintf(The  density  is %1.4f  kg/mA3W,rho2) 

fprintf(The  temperature  is %3.3g  K\n',T2) 

fprintf(The  speed  of  sound  is %3.4f  m/s\n',c2) 

£?  •    l "/  I^^J*"J^"J'»J''J'«i'iJ'*l'»J'«J'  •I*  »1-  *!»  «4»  «A«  -t"  »1*  •1'  •!•  •!•  •!*  -1*  -A-  «1»  »1»  tl>  »1*  •!•  »L»  •!*  •!•  «1»  »1*  •!•  »!•  •!•  »I*  -t»  •!•  •!•  -I*  •!•  •!•  -X»  -1*  »t»  -I-  •!»  - 1»  -l-  -t-  -l*  •!*  *1*  >-U  «1*  vU  •.!*  vL.  ^L.  .L.  k(,  .  L. 

T-|"\-|*|  »"\  tTI       T*  f  T*  *T*  1*  *P  *T*  *T*  1*  *T*  *I*  *T*  *T»  T*  *T*  *T*  *T*  *T*  "T*  *T*  *P  •(*  T*  *P  *?■  *f'  *T*  *T*  *T*  •f  *T*  *T*  I*  *T*  *T*  ^*  *P  *l*  "1*  *T*  *!*  *I*  ^*  1*  *^  *P  T*  T*  T*  T*  *1*  T*  ^*  <P  1*  *P  *P  *|*  *T*  I*  *T*  ■!■  f*  I*  *P 

*W) 

% 

% 

%                                                                      Define  initial  conditions 
% 

for  x=l:stop 

for  y=l:stop 

op(x,y)=Pin; 

01r(x,y)=rho; 

Olru(x,y)=0; 
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end 
end 


Olrv(x,y)=0; 
oE(x,y)=Pin/gl; 


% 

%                      Define  boundary  conditions  at  grid  point  x=l  y=l 
% 

y=2; 

x=2; 

op(x,y)=Ps; 

01r(x,y)=rho/(  1  -(2/(g+ 1  ))*(  1  -( 1  /MachA2))); 
01ru(x,y)=01r(x,y)*Vin; 
01rv(x,y)=01r(x,y)*Vin; 

oE(x,y)=(Ps/g  1  )+(0. 5*(01m(x,y)A2+01rv(x,y)A2)/. 
01r(x,y)); 


% 

%  Calculate  the  new  speed  of  sound  based  on  new  pressure 

%  and  density  values 

% 

nc=sqrt(g*op(x,y)/01r(x,y)); 
% 

%                                                                      Calculate  the  initial  time  step 
% 

dt=CFL*(  1  /(abs(01ru(  1 , 1  )/01r(  1 , 1  ))/dx+abs(01rv(  1 , 1  )/01r(  1,1))/... 
dy+nc*(  1  /dxA2 + 1  /dyA2)  A(0. 5))); 


fprintf(The  inital  time  step  is %1.5es\n',dt) 

% 

%                                                                      Initialize  total  time  counter 
% 

tt=0; 
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% 

%                      Finite  difference  for  density,  momentum  and  energy 
% 

% 

%                     Adjust  x  and  y  grid  end  point  for  finite  differencing 
% 

stop=stop- 1 ; 


for  t=  1 :  endt 

% 

%  NBS  for  column  1  and  row  1 

% 

dtmin=l.; 

dtdx=dt/dx; 
dtdy=dt/dy; 

% 

%                                                                      Calculate  total  time 
% 

tt=tt+dt; 

% 

for  x=l:  stop 
for  y=l:  stop 

% 

%  PREDICTOR 

% 

%  Solve  as  a  function  of  time  and  position  for 

%  increasing  time 
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%- 

% 

%- 


Density(r) 


%- 

% 

%- 


br(x,y)=01r(x,y)-... 
(dt*(((01ru(x+l,y)-01ru(x,y))/dx)+... 

((01rv(x,y+l)-01rv(x,y))/dy))); 

br(x,stop+ 1  )=br(x,stop); 
br(stop+ 1  ,y)=br(stop,y); 
br(stop+ 1  ,stop+ 1  )=br(stop,stop); 


Momentum(ru) 


%- 

% 
%- 


bru(x,y)=01ru(x,y  )-dt* . . . 
((1/dx*... 

(01ru(x+ 1  ,y)A2/01r(x+ 1  ,y)-. . . 
01ru(x,y)A2  /01r(x,y)  +... 
(op(x+l,y)-op(x,y)  '      )))+... 
(1/dy*... 

(01rv(x,y+ 1  )*01ru(x,y+ 1  )/01r(x,y+ 1  )- 
01rv(x,y)  *01ru(x,y)  /01r(x,y)    ))); 

bru(x,stop+ 1  )=bru(x,stop); 
bru(stop+ 1  ,y  )=bru(stop,y); 
bru(stop+ 1  ,stop+ 1  )=bru(stop,stop); 


Momentum(rv) 


brv(x,y)=01rv(x,y)-dt*. . . 
((1/dx*... 

(01ru(x+ 1  ,y)*01rv(x+ 1  ,y)/01r(x+ 1  ,y)- 
01ru(x,y)  *01rv(x,y)  /01r(x,y)  ))+... 
(1/dy*... 

(01rv(x,y+ 1  )A2/01r(x,y+ 1 )-. . . 
01rv(x,y)A2  /01r(x,y)+... 
(op(x,y+l)-op(x,y)       )))); 

brv(stop+l  ,y)=brv(stop,y); 
brv(x,stop+ 1  )=brv(x,stop); 
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brv(stop+ 1  ,stop+ 1  )=brv(stop,stop); 
% 

%  Energy(e) 

% 

bE(x,y)=oE(x,y)-dt*... 
((1/dx*... 

((oE(x+ 1  ,y)+op(x+ 1  ,y))*01ru(x+ 1  ,y )/01r(x+ 1  ,y)  -. . . 
(oE(x,y)+op(x,y)     )*01ru(x,y)  /01r(x,y)    )+... 
1/dy*... 
((oE(x,y+l)+op(x,y+l))*01rv(x,y+l)/01r(x,y+l)  -... 
(oE(x,y)+op(x,y)     )*01rv(x,y)  /01r(x,y)    ))), 

bE(stop+l  ,y)=bE(stop,y); 
bE(x,stop+ 1  )=bE(x,stop); 
bE(stop+ 1  ,stop+ 1  )=bE(stop,stop); 
% 

%  return  for  next  x 

end 

%  return  for  next  y 

end 

% 

for  x=2:stop+l 

for  y=2:stop+l 

% 

%  CORRECTOR 

% 

%  Density(r) 

% 

nr(x,y)=(br(x,y)+01r(x,y))/2-... 
(  dtdx/2*(bru(x,y)-bru(x-l,y))+... 
dtdy/2*(brv(x,y)-brv(x,y- 1 ))); 


nr(l,y)=nr(2,y); 
nr(x,l)=nr(x,2); 
nr(2,2)=01r(2,2); 
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nr(l,l)=01r(2,2); 
%  Left  wall  boundary  conditions 


%- 


%  Momentum(ru) 

% 

Px=gl*(bE(x,y)-0.5*((bru(x,y)A2/br(x,y))+(brv(x,y)A2/br(x,y)))); 

Pxm=g  1  *(bE(x- 1  ,y  )-0. 5  *((bru(x- 1  ,y)A2/br(x- 1  ,y))+(brv(x- 1  ,y)A2/br(x- 1  ,y)))); 

nru(x,y)=(bru(x,y)+01ru(x,y))/2-... 
((dtdx/2*... 

((bru(x,y)  *bru(x,y)  /br(x,y)  +Px)-... 
( bru(x- 1  ,y)*bru(x- 1  ,y)/br(x- 1  ,y)+Pxm)))+. . . 
(dtdy/2*... 

(  bru(x,y)  *brv(x,y)  /br(x,y)-... 
bru(x,y-l)*brv(x,y-l)/br(x,y-l)      ))) ; 

% 


nru(l,y)=nru(2,y); 
nru(x,l)=nru(x,2); 
nru(2,2)=01ru(2,2); 
nru(l,l)=01ru(2,2); 

%  Momentum(rv) 

% 

Py=gl*(bE(x,y)-0.5*((bru(x,y)A2/br(x,y))+(brv(x,y)A2/br(x,y)))); 

Pym=g  1  *(bE(x,y- 1  )-0. 5  *((bru(x,y- 1  )A2/br(x,y- 1  ))+(brv(x,y- 1 )  A2/br(x,y- 1 )))); 

nrv(x,y)=(brv(x,y)+01rv(x,y))/2-... 
((dtdx/2*... 
(  bru(x,y)  *brv(x,y)  /br(x,y) 

bru(x- 1  ,y )*brv(x- 1  ,y )/br(x- 1  ,y) ))+... 
(dtdy/2*... 
( (brv(x,y)*brv(x,y)/br(x,y)+      Py  )-... 

(brv(x,y- 1 )  *brv(x,y- 1  )/br(x,y- 1  )+Pym)))); 

% 
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nrv(l,y)=nrv(2,y); 
nrv(x,l)=nrv(x,2); 
nrv(2,2)=01rv(2,2); 
nrv(l,l)=01rv(2,2); 

%  Left  wall  boundary  conditions 

% 

%  Energy  (E) 

% 

nE(x,y)=(bE(x,y)+oE(x,y))/2-.. . 

((dtdx/2*... 

((bru(x,y)  /br(x,y)  *(bE(x,y)  +Px  ))    -... 

(bru(x- 1  ,y)/br(x- 1  ,y)*(bE(x- 1  ,y)+Pxm))))  +. . . 

(dtdy/2*... 

((brv(x,y)  /br(x,y)  *(bE(x,y)+Py    ))    -... 

(brv(x,y- 1  )/br(x,y- 1 ) *(bE(x,y- 1  )+Pym))))    ); 
% 
% 


nE(l,y)=nE(2,y); 

nE(x,l)=nE(x,2); 

nE(2,2)=oE(2,2); 

nE(l,l)=oE(2,2); 
% 

%return  for  next  x 
end 
% 

%return  for  next  y 

end 
% 

for  y=l:stop+l 
for  x=l:stop+l 
% 

%                      Calculate  pressure  using  the  continuity  equation  for 
%                       the  ideal  gas  law  and  the  new  values  for  energy, 
%                       momentum  and  density 
% 
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np(x,y)=(gl)*(nE(x,y)-(.5*((nru(x,y)A2/. 
nr(x,y))+(nrv(x,y)A2/nr(x,y))))); 


% 

%  Calculate  velocity  using  the  new  values  for  momentum 

%  and  density 

% 

nu(x,y)=nru(x,y)/nr(x,y); 
nv(x,y)=nrv(x,y)/nr(x,y); 


% 

%  Check  stability 

%                                                                      Calculate  the  new  speed  of  sound 
% 

nc=sqrt(g*np(x,y)/nr(x,y)); 
% 

% 

%  Calculate  the  new  time  step  dt 

% 


dtn=CFL*(  1  /(abs(nu(x,y))/dx+abs(nv(x,y  ))/dy+. 

nc*(  1  /dxA2+ 1  /dyA2)A(0. 5))); 
if  dtn  <  dtmin 

dtmin=dtn; 
%  end  dtmin  if  statement 

end 
% 

% 

% 

%retum  for  next  x 

end 
% 

%return  for  next  y 
end 


86 
% 

% 

%                     Update  time  step  and  CFL(Courant  number) 
% 

dt=dtmin; 
fprintf('delta  t=  %1.5e  s\tTotal  time-  %1.5e  s\n',dt,tt) 


% 

%  Output  the  values  for  pressure,  density  and  velocity 

%  using  graphs 

% 

%  Pressure 

for  x=2:stop+l 

for  y=2:stop+l 

press(x-l,y-l)=np(x,y); 
end 
end 

colormap(cool) 
surf  (press) 
hold  on 

pcolor((press)-Pin) 
%contour(press,20) 
title('Absolute  Pressure') 
xlabel('X  grid  points') 
ylabel('Y  grid  points') 
zlabel('Pascals') 

pause(.Ol) 
%hold  off 
% 

%  Replace  old  values  with  new  values  and  begin 

%  next  time  step 

% 

op=np; 


87 


01r=nr; 

01ru=nru; 

01rv=nrv, 

oE=nE; 

% 

%  Return  for  next  t 

end 

% 
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